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Chapter  I 
Introduction 

The  electrical  engineer  must  often  solve  linear,  constant-coefficient,  differential 
equations,  especially  in  such  fields  as  circuit  analysis  and  control  theory.  The  usual 
method  of  solution  involves  use  of  the  Laplace  transform.  Although  this  method 
of  solution  is  well-known  [1],  computing  solutions  to  such  equations  often  involves 
many  tedious  calculations.  The  most  difficult  step  in  this  method  of  solution  is 
performing  the  inverse  Laplace  transform  of  a  rational  function.  The  object  of  this 
thesis  is  to  describe  an  algorithm  for  solving  large  problems  of  this  kind. 

Solving  such  equations  is  not  necessarily  a  numerical  analysis  problem.  We  are 
not  trying  to  approximate  the  solution  to  a  differential  equation.  We  already  know 
the  general  form  of  the  solution,  and  are  simply  seeking  coefficients  for  the  solution. 
Thus  the  problem  is  actually  algebraic,  rather  than  analytic. 

There  are  three  distinct  steps  in  computing  the  inverse  Laplace  transform  of  a 
rational  function.  They  are:  factoring  the  denominator  polynomial,  finding  the  par- 
tial fraction  expansion  of  the  rational  function,  and  computing  the  inverse  Laplace 
transform  of  each  of  these  partial  fractions.  Except  for  the  factoring,  these  are  not 
analytic  problems,  and  even  the  factoring  algorithm  presented  here  makes  use  of  an 
algebraic  technique  to  speed  up  the  finding  of  multiple  roots. 

Special  concern  was  devoted  to  the  partial  fraction  expansion  portion  of  the 
algorithm  in  order  to  reduce  the  number  of  calculations.  The  method  presented 
herein  is  based  on  one  developed  elsewhere  [4]  and  modified  to  further  reduce  the 
number  of  calculations.  The  effort  to  minimize  the  number  of  calculations  required 
constitutes  the  main  contribution  of  this  research. 

As  the  title  declares,  this  thesis  concerns  itself  with  computing  the  inverse 
Laplace  transform  of  a  rational  function,  but  not  with  obtaining  the  rational  func- 
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tion  in  the  first  place.  The  algorithm  consists  of  four  parts.  The  first  part  of  the 
algorithm  accepts  a  rational  function  as  input  and  outputs  the  rational  function 
with  the  denominator  expressed  as  a  product  of  irreducible  factors.  The  next  part 
accepts  as  input  a  rational  function  with  a  factored  denominator.  It  outputs  the 
partial  fraction  expansion  of  the  given  rational  function.  The  third  stage  accepts 
this  expansion  and  computes  the  inverse  Laplace  transform.  The  final  stage  evalu- 
ates the  inverse  transform  over  a  desired  interval  and  plots  a  graph,  if  desired. 


Chapter  II 
Some  Preliminaries 

We  first  review  the  definition  of  the  Laplace  transform  and  its  usefulness  in 
solving  linear,  const  ant -coefficient,  differential  equations.  Given  a  real-valued  func- 
tion /(<)  of  a  real  variable,  its  one-sided  Laplace  transform,  F(s),  is  given  by  [1]: 

•H-oo 


/•-l-oo 

F(s)  =  £{f(t)}  =   /         f(t)e-stdt, 
Jo 


where  5  is  complex  and  hence  F(s)  is,  in  general  complex- valued.  This  transform 
has  two  important  properties  which  can  be  used  to  transform  a  differential  equation 
in  t  into  an  algebraic  equation  in  s.  Both  of  these  properties  are  easy  to  derive  and 
proofs  are  given  elsewhere.  These  properties  are  [1]: 

Linearity.  Suppose  fi(t)  and  f2{t)  are  such  that  their  Laplace  transforms,  Fi(s) 
and  F2(s),  exist.  Also,  let  ci,C2  E  R»  Then 

£{cifi(t)  +  c2/2(r)}  =  c.F^s)  +  c2F2(s). 

Forward  Derivative  Property.  Suppose  that  /(r)  is  (n  —  1)  times  continuously 
differentiable  and  that  the  nth  derivative,  f^n'(t),  is  such  that  its  Laplace  transform 
exists.  Then 

n 

£{/<">(*)}  =  snF(s)  -  53*"~i/°""1)(0). 
Given  an  nth-order  linear  differential  equation  with  constant  coefficients 

n 

]Ta;/(i)(*)  =  K*), 

and  the  initial  conditions,  /(0)  =  c0,/(1^(0)  =  C\, . . . ,  /(n-1)(0)  =  c„_i,  we  can 
take  the  Laplace  transform  of  both  sides  of  the  equation  and  end  up  with  [1] 

F(s)  =         }      i*A        ^    3 ,  (2.1) 

Lj=0  ajsJ 
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where  B(s)  =  C{b(t)}. 

We  see  that  if  b(t)  is  such  that  B(s)  is  a  rational  function,  then  F(s)  will  also 
be  a  rational  function.  Thus,  the  problem  is  reduced  to  finding  the  inverse  Laplace 
transform  of  F(s).  What  might  be  done  when  b(t)  is  such  that  B(s)  is  not  a  rational 
function  will  be  discussed  later. 

The  method  employed  here  to  find  the  inverse  Laplace  transform  of  (2.1)  is, 
first,  to  factor  the  denominator;  next,  to  expand  the  rational  function  into  partial 
fractions;  and  finally,  by  applying  the  linearity  property  of  the  inverse  Laplace 
transform,  to  find  the  inverse  Laplace  transform  of  each  of  the  partial  fractions.  The 
algorithms  presented  in  the  next  three  chapters  accomplish  each  of  these  steps  by 
using  elementary  concepts  of  polynomial  ring  theory,  and  recursion  where  possible. 


Chapter  III 
Factoring 

The  form  of  the  inverse  Laplace  transform  of  a  rational  function  depends  on 
the  factors  in  the  denominator  polynomial.  The  problem  of  factoring  polynomials 
is  an  important  area  of  numerical  analysis.  The  algorithm  presented  here  relies  on 
elementary  algebraic  considerations.  It  should  be  pointed  out  that  any  other  fac- 
toring algorithm  could  be  used  with  the  rest  of  the  programs  listed  in  the  appendix 
by  providing  suitable  interfacing  software. 

The  algorithm  presented  here  is  based  on  the  root-finding  method  of  D.  E. 
Muller  [2],  a  brief  description  of  which  follows.  Given  a  polynomial,  p,  start  with 
three  initial  estimates  for  a  root:  x_2,  x_i,  and  xo-  At  the  nth  stage  (0  <  n  £  N),  fit 
a  quadratic  equation  to  the  points  (xn_2,p(xn_2)),  (in-hP(2n-i))i  and  (xn,p(xn)). 
Then  find  the  root  of  this  quadratic  equation  nearest  xn.  This  root  becomes  xn+i. 
Repeat  this  procedure  while  n  is  less  than  a  given  upper  bound,  or  subsequent 
estimates  fail  to  improve  by  a  given  small  amount.  To  find  other  roots,  deflate  p  by 
the  appropriate  factor  and  reapply  Muller's  method  if  necessary. 

Suppose  we  are  given  a  polynomial  with  real  coefficients  and  distinct  complex 
roots.  Apply  Muller's  method  to  obtain  a  root.  Deflate  the  polynomial  by  a  factor 
containing  the  root  just  found  (and  its  complex  conjugate,  if  necessary)  to  obtain 
a  polynomial  of  smaller  degree  with  real  coefficients  and  distinct  complex  roots. 
Repeat  this  process  until  the  degree  of  the  deflated  polynomial  is  of  degree  0,  at 
which  point  all  the  roots  of  the  original  polynomial  have  been  found. 

The  above  process  works  quite  well  if  all  the  roots  are  distinct.  If  Muller's 
method  is  applied  to  find  a  multiple  root,  it  converges  more  more  slowly  than  it 
does  for  a  single  root.  That  is,  it  runs  through  many  more  iterations  before  the 
difference  between  successive  estimates  becomes  as  small.    So  we  turn  to  abstract 


algebra  for  a  way  to  speed  things  up.  By  using  the  method  described  below,  we  can 
ensure  that  we  will  always  be  searching  for  distinct  roots. 

Consider  a  monic  polynomial  p  E  R[x],  with  r  £  C  a  root  of  p.  Then  there 
exists  q  £  C[x]  such  that  p  =  (x  —  r)q.  Now  consider  the  first  derivative  of  p,  p'  = 
q  +  (x  —  r)q' .  Observe  that  p'{r)  =  0  if  and  only  if  q(r)  =  0.  Thus  r  appears  more 
than  once  as  a  root  of  p  if  and  only  if  p'(r)  =  0  [3]. 

Now  look  at  gi  =  gcd(p,p')  €  R[x].  If  deg(gi)  =  0  then  the  roots  of  p  are  all 
distinct.  If  deg(^i)  >  0  then  each  root  of  g\  is  also  a  root  of  p  and  appears  more 
than  once  as  a  root  of  p.  Now  define  gn  =  gcd(p,//n)).  If  deg(gn)  >  0  then  each 
root  of  gn  is  also  a  root  of  p  and  its  multiplicity  is  at  least  n  +  1. 

Since  the  multiplicity  of  any  root  is  at  most  deg(p)  there  exists  a  minimal  k  €  N 
such  that  deg(<7fc)  =  0.  This  k  can  be  found  by  repeated  differentiation  of  p  and 
application  of  the  Euclidean  algorithm  to  obtain  each  gcd.  The  first  gcd  obtained 
with  degree  0  is  <?*. 

The  only  roots  of  gk-i  will  be  all  the  roots  of  p  with  multiplicity  k.  These 
roots  can  be  found  using  Muller's  method  on  gk-i-  These  roots  will  be  distinct 
in  gk-i  but  of  multiplicity  two  in  gk-2-  But  the  remaining  roots  of  gk-2  will  be 
distinct,  and  can  be  found  by  deflating  gk-2  by  all  known  roots  and  then  applying 
Muller's  method  on  the  deflated  polynomial.  Similarly,  all  roots  of  p  of  multiplicity 
n  or  higher  will  be  roots  of  gn-\-  Say  a  root  has  multiplicity  m  >  n;  this  root  will 
have  multiplicity  m  —  n  +  1  in  gn-\.  Thus,  we  can  find  all  roots  of  p  by  working 
backward  from  gk  to  p  using  Muller's  method  and  deflation.  Most  importantly, 
Muller's  method  is  never  used  to  search  for  a  root  of  multiplicity  greater  than  1. 

It  should  be  pointed  out  that  repeated  polynomial  divisions  occur  in  the  Eu- 
clidean algorithm.  A  loss  of  precision  is  inherent  in  this  process.  If  one  cannot  make 
this  sacrifice  of  accuracy  of  the  results  in  favor  of  increased  rate  of  convergence,  re- 
alize that  any  other  factoring  algorithm  will  work  with  the  rest  of  the  programs  in 
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the  appendix.  The  interfacing  software  would  be  easy  to  write. 


Chapter  IV 
Partial  Fraction  Expansion 

One  objective  in  designing  a  partial  fraction  expansion  algorithm  was  to  min- 
imize the  amount  of  calculation  done.  Chin  and  Steiglitz  [4]  devised  an  algorithm 
capable  of  accomplishing  the  expansion  in  N(N  —  1)  multiplications  and  |iV(  JV  —  1) 
additions,  where  N  is  the  degree  of  the  denominator  of  the  given  rational  function. 
This  algorithm  has  a  disadvantage:  it  requires  use  of  complex  arithmetic.  Chin  and 
Steiglitz  count  complex  divisions  as  equivalent  in  time  to  complex  multiplications. 
While  this  may  be  true  for  real  divisions  and  multiplications  it  certainly  is  not  true 
for  complex  ones.  Observe  that, 

a  +  ib       (ac  +  bd)  +  i(bc  —  ad) 
c  +  id  '  c2  +d2  ' 

has  6  real  multiplications,  2  real  divisions  and  3  real  additions.   Call  this  8  multi- 
plications and  3  additions.  Furthermore, 

(a  +  ib)(c  +  id)  =  (ac  —  bd)  +  i(bc  +  ad), 

requires  4  real  multiplications  and  2  real  additions.  Finally,  note  that 

(a  +  ib)  +  (c  +  id)  =  (a  +  c)  +  i(b  +  d), 

consists  of  2  real  additions. 

Examination  of  Chin's  and  Steiglitz's  algorithm  reveals  that  the  expansion 
actually  involves  |iV(iV-l)  complex  additions,  |JV(JV—  1)  complex  multiplications, 
and  ^N(N  —  1)  complex  divisions.  Using  the  above  calculations,  this  results  in 
6N(N  -  1)  real  multiplications  and  ^A^iV  -  1)  additions.  So  if  complex  arithmetic 
can  be  avoided  and  the  operation  count  can  be  held  lower  than  this,  the  algorithm 
will  be  improved. 
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It  is  clear  that  the  reason  Chin  and  Steiglitz  chose  to  work  with  complex  num- 
bers is  that  a  polynomial  in  R[x]  C  C[x],  splits  in  C.  This  makes  the  partial  fraction 
expansion  algorithm  simple  to  describe  and  analyze.  But  R[x]  has  the  property  that 
an  irreducible  element  is  either  linear  or  quadratic.  If  we  use  this  property  we  can 
avoid  complex  arithmetic  and  thus  reduce  the  number  of  calculations  at  the  cost  of 
complicating  the  algorithm  a  bit. 

The  key  to  adapting  the  algorithm  of  Chin  and  Steiglitz  is  finding  a  nice  way 
to  generalize  the  following  problem.  Let  s,  r,  K  £  C,  for  s/r  and  find  A,  A*  £  C, 
such  that 


x  —  s  \_(x  —  r)n 


A* 


+ 


A 


(x  -  s)(x  -  r)n_1        (x-r)n 
This  problem  generalizes  to:    Let  s,r  £  R[x]  \  R,  with  gcd(.s,r)  =  1;  and  A'  £ 
R[x],  with  deg(k)  <  deg(r).    Find  A,  A*  £  R[x]  such  that  deg(A)  <  deg(r),  and 

deg(.4*)  <  deg(.s),  and 

1  T  K 1  A  *  A 

(4.1) 


r  +  — 

srn-\  rn 


First  multiply  through  by  rn       to  reduce  the  problem  to: 


K 

r 


A*       A 

s  r 


(4.2) 


We  know  that  since  gcd(.s,r)  =  1,  there  exists  a,  b  £  R[x]  such  that  as  +  br  =  1 

and  so  K  =  K[as  +  br]  =  Kas  +  Kbr.  Apply  the  division  algorithm  to  obtain  the 

following: 

Ka  =  rq  +  A     such  that      deg(.4)  <  deg(r), 

(4.3) 


Ka  =  rq  +  A     such  that      deg(.4)  <  deg(r), 
Kb  =  sq*  +  A"     such  that      deg(.4*)  <  deg(^). 


Now  write 


K  =  Kas  +  Kbr  =  (rq  +  A)s  +  (sq*  +  A*)r 

=  (q  +  q*)rs  +  As  +  Amr. 
By  hypothesis  we  have  deg(if)  <  deg(r)  <  deg(r.s)  and  from  (4.3)  we  get  deg(.4.s)  < 
deg(r.s)  and  deg(.4*r)  <  deg(rs).   So  q  -f  q*  =0  and  K  =  As  +  Amr,  which  yields 
(4.2)  as  required. 


Now  we  must  determine  how  to  calculate  A  and  A*  in  (4.1)  [5],  and  the  neces- 
sary number  of  calculations.  First  consider  the  following  adaption  of  the  Euclidean 
algorithm.  Let  s,r  £  R[x].  The  division  algorithm  gives: 

s  =  rq1+  xj  deg(xi)  <  deg(r) 

r  =  xiq2  +  x2  deg(x2)  <  deg(xi) 

xi  =  Z2?3  +  2?3  deg(x3)  <  deg(x2) 

xn_2  =  xn_j5n  +  xn  deg(xn)  <  deg(xn_i) 
And  also  define  xq  =  r.  If,  furthermore: 

a0  =  0  b0  =  1 

ai  =  I  &i  =  -qi 

an  =  O-n-2  —  qnO-n-l  bn   =  &n-2  ~  <Zn&n-l 

then  an3  +  bnr  =  xn,  for  all  n  >  0. 

Proof:  ao-s  +  &or  =  r  =  xo-  a\S  +  b\r  =  s  —  rqi  =  x\.  Assume  the  hypothesis 
is  true  for  n  —  2  and  n  —  1  for  some  n  G  N. 

^n  =  ^n-2  ~  Qnxn  —  1 

=  an-2>s  +  6n-2^  -  ?n(an-i-s  -I-  bn-ir) 
=  (an_2  -  qnan-i)s  +  (6n-2  -  9n&n-i)r 
=  ans  +  6nr, 

as  claimed. 

We  can  use  the  above  procedure  to  develop  a  method  for  evaluating  A  and 
—A*.  There  are  four  cases  to  consider  since  either  s  or  r  can  be  of  degree  1  or  2. 
Let  us  begin  with  the  most  difficult  case: 

Case  1:   Both  5  and  r  are  quadratic 
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We  have  s  =  rqi  +  x  and  r  =  X\q2  -f  x2.  x2  is  a  unit  so  a2s  +  b2r  =  —q2s  -f 
(1  +  q\q2)r  =  x2,  and  we  can  get  the  following. 

K        ±K(-q2s  +  (l  +  qiq2)r) 
sr  sr 

r  5 

By  the  division  algorithm,  92 -K"  =  Qr  —  x2A  for  some  Q  G  R[x].  The  quotient, 
Q,  does  not  matter  as  was  shown  in  the  derivation  of  (4.1).  The  remainder  is  the 
important  thing.  Now  we  need  to  add  up  all  the  calculations  required  to  evaluate 
A. 

Table  (4.1)     Summary  of  Case  1 

To  obtain  x\  requires  2  adds, 

to  obtain  q2  and  x2  requires  2  divs     2  mult  2  adds, 

to  obtain  q2K  requires  4  mults  2  adds, 

to  divide  by  r  to  get  x2A  requires  2  mults  2  adds, 

and  to  evaluate  A  requres  2  divs     1  add. 

Resulting  in  4  divs     8  mults  9  adds. 

Now  we  have  A.  We  need  —A*  as  well.  Since  we  know  that  K  =  As  +  A*r,  write 
K  =  kix  +  ko,  A  =  ajx+ao?  A*  =  a*x  +  cto>  s  —  x2  +SiX  +  s0,  and  r  —  x2  +rix  +  r0. 
We  obtain 

k\x  +  k0  =  (a\x  +  clq)(x2  +  8\X  +  So)  +  (a*x  +  al){x~  +  ri&  +  ro). 

Equating  coefficients  for  cubes  and  constants  results  in  —a*  =  at ,  and  — Oq  = 
(ciqSo  —  ko)/r0.  So  —.4*  can  be  calculated  using  one  multiplication,  one  division, 
and  one  addition.  All  told,  calculation  of  .4  and  —A*  requires  the  equivalent  of  14 
multiplications  and  11  additions. 
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Case  2:  s  is  quadratic  and  r  is  linear 

We  have  5  =  rqx  +  x\.  Since  x1  is  a  unit,  we  stop. 

K       j^K{s  -  qir) 


sr  sr 


r  s 

—  =  A,  and  both  K  and  xi  are  units,  so 

Table  (4.2)     Summary  of  Case  2 

to  obtain  X\  requires  1  mult       2  adds, 

and  to  obtain  A  requires     1  div.     

The  result  is  1  div      1  mult       2  adds. 


Now  we  have,  as  before, 

kxx  +  kQ  =  A(x2  +  s\x  +  so)  +  {a\x  +  a%){x  +  r0). 

Equating  coefficients  of  squares  and  constants  gives  —a*  =  A  and  —  clq  =  (Asq  — 
ko)/ro-  So  —  A*  is  computed  with  one  multiplication,  one  division,  and  one  addition. 
All  together,  calculation  of  A  and  —A*  requires  the  equivalent  of  4  multiplications 
and  3  additions. 

Case  3:  s  is  linear  and  r  is  quadratic 

We  have  s  =  rq\  +  xi,  and  r  =  x\q2  +  X2-  But  q\  =  0,  thus  x\  =  s  and  notice 
that  X2  is  then  a  unit,  so  —  ^s  +  r  =  £2  and  then 

A"  _  j;K(-q2s  +  r) 
sr  sr 

_  X2^Z  I       *2 


r  s 

This  time  we  will  compute  —A*,   a  unit,  first.     Observe  that  K   =   Qs  — 
(x2(—A*)).     Now  we  do  as  before  and  write  k\x  +  k0   =   [a\X  +  a0)(x  +  s0)  + 
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A*(x2  +  rxx  +  r0).    Equate  coefficients  of  squares  and  constants  to  get  ax  =  —A* 
and  a0  =  (k0  —  A*tq)/sq.  Now  we  summarize. 

Table  (4.3)     Summary  of  Case  3 

To  obtain  x2  requires  1  mult  2  adds, 

to  obtain  x2 A*  requires  1  mult  2  adds, 

to  obtain  —  A*  requires  1  div  1  add, 

and  to  obtain  a0  requires  1  div       1  mult  1  add. 

The  result  is  2  divs     3  mults  5  adds. 

All  together  calculation  of  A  and  —A*  requires  the  equivalent  of  5  multipli- 
cations and  5  additions. 

Case  4:  r  and  s  are  both  linear 

This  case  is,  of  course,  the  simplest,  s  =  q\r  +  X\,  x\  is  a  unit  and  q\  =  1.  So 

K       ±K{a-r) 

sr  sr 

K       —K 

-  JlL  _i_        Xl 


r  s 

Thus  A  =  —A*  =  iv/xi,  which  requires  1  multiplication  and  1  addition.  Now 
that  each  case  has  been  examined,  the  following  table  summarizes  the  preceding 
information: 
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Table  (4.4)     Number  of  operations  required  to  evaluate  (4.2) 

deg(s)     deg(r)     multiplications  additions 
111  1 

2  1  4  3 

12  5  5 

2  2  14  11 

The  following  is  essentially  Chin's  and  Steiglitz's  algorithm  in  R[x]  instead  of 
C[x].  Let  p  €  N  be  given  and  dj  G  R[x],  for  each  1  <  j  <  p  and  each  dj  irreducible 
over  R.  Let  rrij  G  N  denote  the  multiplicity  of  each  dj.  Also  let  Qo,D  £  TL[x]  be 
given  such  that  D  =  YVj=i  (^j)m;  and  deg(Q0)  <  deg(Z?)  =  5Zi=i  deg(<ij)m;.  Thus 
we  have  a  proper  rational  function  and  wish  to  find  K{j  £  R[x]  such  that 

Define  mo  =  0  and  n  =  ^2Pj-0  rrij.  Now,  for  1  <  /  <  n,  define 

u 


1  =  1  + 

min{x  G  N  U 

{0}lX>>< 

j=0 

a, 

•J- 

if  j  <  u/; 
rrn,     if  j  =  u/, 

for  I  <  j   <  uj.    Also,  for  each  /,  define  //  =  dU/,  and  Ri,Qi  G  R[x]  such  that 

Qj-i  =  Q///  +  Ri.  Again  define  A0(/v),  A^(K )  :  R[x]  h*  K[x]  by 

tf    _  Af/A')       A0(7v) 


Now  the  partial  fraction  expansion  can  be  obtained  as  below  in  n  steps,  the 
Ith  step  being: 

Q°  1     -IQo] 


nj_, <"    n-=1  a 


A  A  A"  (4-5) 


■K'  +  EEtt^ 
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where  it  is  understood  that  Iln>Jfc>n  fk  =  ^  and 


[0, 

4-     W'  +  ^wy)). 


K 


1-1 
(i-i)i' 


if  i  >  uj  or  j  >  m  or  I  =  0; 
if  j  <  u/  and  i  =  v1.; 
if  j  <  m  and  i  <  uj; 
if  j  =  if/  and  i  >  1; 


/-i 


I  Ri  +  £o<*<i  ^  W;1  +  AikWk)),     *  J  =  «i  and  i  =  1 


(4.6) 


Now  we  must  count  operations  to  compare  this  method  with  Chin's  and  Stei- 
glitz's  method.  It  turns  out  that  the  number  of  operations  depends  on  the  number 
of  quadratic  factors  in  the  denominator  of  the  rational  function.  Let  N  =  deg(D) 
and  then  denote  the  number  of  quadratic  factors  in  D  as  9.  Thus  N  =  n  +  q. 

First  consider  the  number  of  calculations  necessary  to  compute  {R(  \  1  <  /  < 
n}.  The  /th  stage  involves  dividing  Qi-\  by  //.  Two  facts  are  necessary:  To  divide 
an  A/-degree  polynomial  by  a  monic  linear  factor  requires  M  multiplications  and  M 
additions  and  to  divide  the  same  polynomial  by  a  monic  quadratic  factor  requires 
2(M  —  1)  multiplications  and  additions.  Note  that  the  largest  that  deg(Qo)  can  be  is 
N  —  1 .  We  shall  prove  that  in  this  worst  case  it  requires  no  more  than  |JV(iV  —  1)  —  q 
multiplications  and  additions  to  compute  {Ri  \  1  <  /  <  n}. 

We  shall  use  induction  on  N.  For  N  =  1  we  must  have  n  —  1,  and  q  =  0. 
and  of  course,  deg(Qo)  =  0  thus  ^N(N  —  1)  —  q  —  0,  which  reflects  the  fact  that 
there  is  really  nothing  to  do  in  this  case.  We  will  also  need  to  examine  the  case 
where  N  =  2,  with  n  =  1  and  q  =  1.  We  still  get  ±N(N  -  1)  -  q  =  0.  which 
again  indicates  that  there  is  really  no  partial  fraction  expansion  to  carry  out.  Now 
assume  the  result  for  a  given  N. 

First  assume  that  we  add  a  linear  factor  to  D  and  increase  deg(Qo)  to  (X  -f 
1)  —  1  =  N.  So  divide  Qq  by  this  new  linear  factor  to  get  deg(Qi)  =  N  —  1,  which 
will  require  N  multiplications  and  additions.  Now  apply  the  induction  hypothesis 
to  Q\.  It  will  require  jN(N  —  1)  —  q  multiplications  and  additions  to  obtain  {Ri  | 
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2  <  /  <  n  +  1}.  Adding  up,  we  get 

N+±N(N-l)-q  =  ±(N  +  l)N-q. 

Now  assume  that  we  add  a  quadratic  factor  to  D  and  increase  deg(Q0)  to 
(N  +  2)  —  1  =  N  +  1.  Dividing  Qo  by  the  new  quadratic  factor  requires  2N 
multiplications  and  additions.  We  are  left  with  deg(Q!)  =  N  —  1.  By  induction,  to 
compute  {Ri  \  2  <  I  <  n}  requires  |iV(JV  —  1)  —  q  multiplications  and  additions. 
Summing,  we  get 

2N  +  \N(N  -l)-q  =  i(iV  +  2)(N  +  1)  -  (q  +  1), 

as  required. 

We  must  also  consider  the  necessary  number  of  calculations  required  to  compute 
{K\i  |  1  <  j '  <  u/,  andl  <  i  <  vlA  for  some  1  <  I  <  n.  It  turns  out  that  the  number 
of  operations  needed  to  compute  this  set  depends  on  deg(dUl),  on  mU(,  and  on  the 
number  of  quadratic  factors  preceding  du, 

Notice  that  computation  of  K\u  requires  no  calculation  for  i  >  1.  This  means 
that  the  largest  operation  count  for  the  partial  fraction  expansion  algorithm  occurs 
when  all  the  factors  of  D  are  distinct,  that  is,  when  m,j  =  1  for  all  1  <  j  <  p  =  n. 
Consequently,  uj  =  /,  so  /j  =  d\  and  Vj  =  rrij  =  1  for  all  admissable  /  and  j.  Now 
write  (4.5)  as 


rijk=l  /*  lln>Jfc>//*  j—i     J  J 

And  we  can  also  write  (4.6)  as 

(  0,  if  /  =  0 

K{j  =  \  ^(Klj1),  ifl<j</ 

l^  +  Eo^i^C^ll1).   if^^- 

Calculation  of  {K[,  \  1  <  j  <  1}  for  a  given  /  requires  that  A\j  and  —  A*j  be 
determined  /  —  1  times  along  with  (/  —  1)  deg(//)  additions.  Consider  the  number  of 
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calculations  in  computation  of  Aij  and  A*-.  Table  (4.4)  gives  us  this  information  if 
we  let  s  =  fi  and  r  =  /,-.  Notice  that,  given  //,  it  will  take  the  most  operations  if 
fj  is  quadratic.  Hence,  the  number  of  calculations  in  computing  {A'{  •  |  1  <  j  <  1} 
will  be  largest  if  deg(/jt)  =  2  for  all  k  <  I. 

In  order  to  find  an  upper  bound  on  the  number  of  calculations  in  this  algorithm, 
assume  all  factors  of  D  are  distinct  and  ordered  such  that  all  quadratic  factors 
appear  first.  Then  the  entire  algorithm  would  require 

1  q 

N(N  _  i)  _  q  +  £  I4(fc  -  1)  +  £(5<?  +  (k  -  (q  +  1)) 

=  N(N  -  1)  -  q2  +  ZNq  -  Iq 

multiplications.  This  formula  also  works  for  q  =  0  and  q  =  N/2.  The  result  in  each 
case  is  N(N  —  1)  and  |iV(JV  —  2),  respectively  which  is  easily  verified.  Also  the 
algorithm  will  require 


I/V(.V-l)-g  +  ]T(ll  +  2)(fc-;L) 


fc=l 

n 


+  £((5  +  2)g  +  (l  +  l)(*-(«  +  l))  ^4-S) 

7+1 
=  lN(N-l)-^q2+3Nq-l-±q 

additions.  Again,  for  the  special  cases  q  =  0  and  q  =  N/2,  the  formula  yields 
%N(N  —  1)  and  ^-N(N  —  2)  respectively.  In  fact  it  is  quite  easy  to  show  that  for 
all  N  >  2  and  N/2  <  q  <  0  we  get  the  following: 

6N(N  -  1)  >  -N(N  -  2)  >  7V(iV  -  1)  -  q2  +  3A'g  -  lq. 

This  shows  that  the  greatest  number  of  multiplications  occurs  when  q  =  n,  and  is 
still  less  than  the  number  required  by  Chin's  and  Steiglitz's  algorithm.  Also  observe 

lAN{N  _  i)  >  hN{N  _  2)  >  In(N  -  1)  -  '-f  +  3iNTg  -  ^g. 
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Which  tells  the  same  story  for  additions. 

These  results  show  that  this  adaption  of  Chin's  and  Steiglitz's  algorithm  saves 
calculations.  To  be  fair,  however,  one  must  realize  that  the  output  of  this  adaptation 
is  not  the  same  as  that  of  Chin  and  Steiglitz.  Which  algorithm  is  better  will  depend 
on  the  application.  Partial  fractions  expansions  can  be  useful  in  a  wide  scope  of 
problems  involving  integrals  of  rational  functions  [6]. 

Chin's  and  Steiglitz's  output  differs  from  the  one  presented  here  in  that  C[x]  is 
the  ambient  polynomial  ring  and  each  denominator  in  the  partial  fractions  expansion 
is  thus  linear.  With  the  adapted  algorithm,  R[x]  is  used  and  the  denominators  can 
be  linear  or  quadratic.  It  happens  that  in  computing  inverse  Laplace  transforms, 
either  form  is  acceptable  and  it  is  better  to  have  fewer  operations;  however,  this 
may  not  always  be  so  for  other  applications  of  partial  fraction  expansion. 
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Chapter  V 
Inverse  Laplace  Transform 

The  most  elementary  approach  to  finding  the  inverse  Laplace  transform  of  a 
given  function  is  to  use  a  table  of  transform  pairs.  Indeed,  large  tables  of  transform 
pairs  have  been  prepared.  In  particular,  a  popular  table  [7]  lists  the  following 
transform  pair: 

fc-l/2 


The  functions  Jp  are  known  as  the  Bessel  functions  of  half-integral  order.  They 
have  the  following  recursive  definition  [8]: 

•7-.i/2(0  =  y  3costi 

Ji/2(t)  =  J—sint, 
Jp+1(t)  =  ^-Jp(t)-Jp.l(t). 

Let  us  simplify  things  somewhat  by  defining: 

fc-l/2 


Hk(at)  =  y/Tri—J  Jk-i/2(at)- 


We  thus  obtain  the  recursive  relationship: 

2 
Ho(at)  =  -  cos(at), 

Hi(at)  =  -  sin(ai), 
a 

2k  -I  /t  \2 

Hk+1(at)  =  -—- J7*(at)  -     —      Hk-i 

la*  \-a/ 

Hence  the  Laplace  transform  pair  (5.1)  becomes: 

1 


(52+a2)M  "  iW 
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If  we  use  a  well-known  property  of  the  Laplace  transform,  we  can  derive  another 
important  Laplace  transform  pair.  Use 


^{^w}  =  -*/(') 


to  get 


•-i 


c-1 


-2ks 


(s2+a2) 

s 
(s2  +  a2) 


fe+i 


fe+i 


T(k) 

t 

W) 

t 


2r(fc  +  i) 


Hk(at) 
Hk(at) 
Hk(at). 


Which  is  equivalent  to 


--i 


■Hk-i(at). 


(s2+a2y  2T(k) 

This  would  be  true  for  k  >  2  but  also  holds  for  k  =  1.  Again  make  use  of  a  Laplace 
transform  property,  namely  C~l  {F(s  +  r)}  =  e~Tt f(t)  to  get  one  of  the  Laplace 
transform  pairs  useful  in  this  problem: 


--1 


((s  +  r)2  +  a2)fc 


■rt 


T(k) 


At 
BHk(at)  +  —Hk-1(at 


(5.2) 


We  also  make  use  of  another  often-tabulated  Laplace  transform  pair  [7]: 

A       1        Ae-ritk-1 


>-i 


(5.3) 


(s  +  r)*J  T(k) 

One  of  these  two  Laplace  transform  pairs  will  apply  to  each  term  in  the  partial 
fraction  expansion  of  a  rational  function.  From  (4.3)  and  we  have 


•  -i 


Qo 


uu  w 


,^-EE 


-  Kij 


tztzw 


j=l   :=1 


-  r-l  J    Ki 


(djY 


[5.4) 
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Now  if  deg(Ktj)  =  1  use  (5.3),  and  if  deg(Kij)  =  2  use  (5.2). 
Using  (5.2),  and  induction  it  can  be  shown  that: 

C'1  I     A(s  +  T)  +  B     I  _  e-rt  y  tj  j      co^at)  +  £,  sin(at)) 
{((s  +  rf+a*)    J  JS 

And  also  note  that  (5.3)  can  be  written 


1  }  — r  }  =  e~Tttk-1  (Acos(O)) 


Thus  we  can  express  (5.4)  in  the  form: 

E  E  £_1     7TT7  \=  E  e_T; '  E  *'  (a*  cos(a^  +  ft*  sin^))  •  <5'5) 

This  is  the  way  the  inverse  Laplace  transform  is  computed  by  the  program  in 
the  appendix.  The  output  is  simply  an  array  of  coefficients  for  an  expression  of  the 
form  (5.5). 
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Chapter  VI 
Applications 

One  of  the  most  probable  applications  of  this  algorithm  will  be  to  evaluate 
transient  responses  of  control  systems  with  a  known  transfer  function.  Figure  (6.1) 
shows  a  control  diagram  for  an  automatic  flight  control  system  for  a  supersonic 
aircraft  [9].  The  transfer  function  corresponding  to  figure  (6.1)  will  depend  on 
the  values  of  K\  and  K-i-  It  is  simple  to  compute  this  transfer  function  using  the 
coefficient  values  shown  in  figure  (6.1). 

The  above  transfer  function  was  inverse  transformed  using  various  values  for 
K\  and  K^.  The  transient  response  functions  so  obtained  are  graphed  in  figures 
(6.2)  and  (6.3).  Notice  how  the  response  improves  until  the  onset  of  instability. 

It  was  mentioned  before  that  sometimes  we  want  the  inverse  transform  of  some- 
thing other  than  a  rational  function.  One  common  example  arises  when  a  control 
system  contains  time  delays.  We  have  the  following  transform  pair. 

e-sT 
£{u(t-T)}  =  


In  order  to  apply  the  algorithm,  we  must  approximate  eaT  by  a  rational  function. 
This  can  be  done  using  the  Pade  approximant  [9].  We  have 


.-.r_.„,.™.    Y.U(-l)'HsT)' 


e-'  ~  Pn(sT)  = 


where 

—  8 

Using  the  Pade  approximant,  the  algorithm  was  used  to  invert  -j-.  Pade 
approximants  of  order  2,3,4,  and  6  were  used.  These  approximants  are  shown 
graphically  in  figure  (6.4). 
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Figure  (6.1)  Control  system  diagram  of  SST  aircraft.  Adapted  (9). 
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Figure  (6.2)  Step  response  of  above  control  system  for  K\Ki  =  0,0.2,0.4. 
and  0.6. 
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Figure  (6.3)  Step  response  of  above  control  system  for  K\ K2  =  0.8,  1.0S,  and 
1.10. 
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Figure  (6.4)  Pade  approximants  to  unit  time  delay. 
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Chapter  VII 
Conclusion 

This  paper  has  presented  an  algorithm  for  computing  inverse  Laplace  trans- 
forms of  rational  functions  as  might  arise  in  practical  electrical  engineering  prob- 
lems. Programs  written  to  demonstrate  the  algorithm  follow  in  the  appendix. 

Numerical  analysis  aspects  of  the  problem  were  not  dealt  with,  but,  except  for 
root-finding,  the  problem  was  shown  to  be  an  algebraic  one.  Results  from  elemen- 
tary abstract  algebra  were  used  to  derive  the  methods  described.  Special  effort  was 
made  to  reduce  the  number  of  calculations  in  the  partial  fraction  expansion. 

Some  applications  were  presented  to  show  practical  results.  These  applications 
also  made  it  clear  that  assuming  that  the  Laplace  domain  function,  F(s),  in  (2.1) 
is  a  ratio  of  polynomials  is  not  always  valid.  Future  work  on  this  problem  should 
concern  itself  with  this  assumption. 
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Appendix  II 
Glossary  of  Terms 

C  denotes  the  field  of  complex  numbers. 

C[x]  denotes  the  ring  of  polynomials  with  coefficients  in  C. 

deg:  If  p(x)  =  do  +  <L\X  +  •  •  •  +  anxn  ^  0  and  an  =fi  0  then  deg(p(x))  =  n. 

Division  algorithm:  Given  two  polynomials  p(x),  q(x)  £  F[x),  where  F[x]  is  the 
ring  of  polynomials  with  coefficients  in  the  field  F  and  q(x)  7^  0,  there  exist  two 
polynomials  t(x),r(x)  £  F[x]  such  that  f(x)  =  t(x)q(x)  +  r(x)  where  r(x)  =  0  or 
deg(r(x))  <  deg(q(x)).  The  process  by  which  t(x)  and  r(x)  are  found  is  known  as 
the  division  algorithm  and  is  simply  the  "long-division"  process  everyone  knows  to 
divide  one  polynomial  by  another. 

Euclidean  algorithm:  Given  two  polynomials  p(x),q(x)  £  F[x],  where  F[x]  is 
the  ring  of  polynomials  with  coefficients  in  the  field  F  and  p(x)  and  q(x)  are  not  both 
0,  then  gcd(p(x),q(x))  E  F[x]  exists  and  there  exist  polynomials  m(x),n(x)  E  F[x] 
such  that  gcd(p(x),q(x))  =  m(x)p(x)  +  n(x)q(x).  The  process  used  to  determine 
these  special  polynomials  is  called  the  Euclidean  algorithm  and  is  shown  explicitly 
in  Chapter  IV. 

T:  A  sufficient  definition  of  the  function  T  for  n€NU{0}is 

JO,  if  n  =  0; 

L[n)~  \(n-l)!,     if  n  >0. 

gcd:  Let  a,  b  £  F[x].  If  c  €  F[x]  satisfies: 

1.  c  is  monic. 

2.  c  divides  a  and  6. 

3.  Any  other  divisor  of  a  and  b  divides  c. 

then  c  is  called  the  greatest  common  divisor  of  a  and  b  and  is  denoted  gcd(a,  b). 

irreducible:  Let  p  £  F[x]  be  such  that  p  =  ab  for  some  a,  6  £  Ffx]  if  and  only 
if  deg(a)  =  0  or  deg(6)  =  0.    Such  a  p  is  said  to  be  irreducible  over  F.  Note  that 
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irreducibility  depends  on  the  field,  F.  x2  +  1  is  irreducible  over  R[x]  but  not  over 
C[x}. 

min  is  a  function  that  operates  on  a  well  ordered  set  and  whose  value  is  the 
minimum  element  of  that  set.  The  well  ordering  property  of  N  asserts  that  if  A  C  N 
then  min  A  exists. 

N  denotes  the  ring  of  natural  numbers. 

R  denotes  the  field  of  real  numbers. 

R[x]  denotes  the  ring  of  polynomials  with  coefficients  in  R.  u(t):  Let  t  €  R. 

fO,     ift<0; 
U(<)-\1,     ift>0. 

unit:  u  6  F[x]  is  called  a  unit  iff  deg(u)=0. 


28 


Appendix  III 
FORTRAN  Programs 

The  following  programs  are  intended  to  merely  demonstrate  the  algorithms 
described  in  the  previous  chapters.  Specifically,  they  were  used  to  evaluate  the 
example  applications  mentioned  in  Chapter  VI.  No  guarantee  of  their  usefulness  to 
any  other  application  is  implied. 

These  programs  could  certainly  be  made  more  user- friendly.  There  are  no  error 
handling  routines,  the  user  interaction  is  minimal,  and  file  management  is  cumber- 
some. Such  things  are  left  to  a  better  programmer.  Nevertheless,  the  programs 
serve  their  purpose  of  demonstrating  the  algorithms. 

The  programs  fall  into  four  slightly  overlapping  categories:    those  associated 

with  Chapters  III  through  V  and  those  for  creating  data  files  and  making  plots. 

The  following  is  a  brief  categorical  index  of  the  programs.  Subroutine  dependence 

is  indicated  by  indentions. 

Factoring  Program 

ROOT_FIND 

POLY_READ 
FACTORER 
DERIV 
EUCLID 

POLDIV 
FIND_EM 

MULLER 

DEFVAL 
COMPOSE 
SPEC-WRITE 
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Partial  Fraction  Expansion  Program 


PARTJFRAC 

SPEC-READ 
EXPAND 

TRANSFER 

POLDIV 

ALIKE 

TRANSFER 
DIFFERENT 

TRANSFER 
EUCLIDEAN 

TRANSFER 
POLDIV 
POLMULT 
POLADD 
PART-WRITE 


Inverse  Laplace  Transform  Program 


INVERT 

INV2 

INIT 

BESSEL 

GAMMA 


PLOT 

READ 

LOTS_0_PLOTS 

PLOT.O_MATIC 


Plotting  Program 


Data  Entry  Routines 

INPUTJIAT 

SPECJNPUT 

SPEC-WRITE 

The  programs  that  follow  are  listed  in  alphabetical  order  according  to  their 

VAX  FORTRAN  filenames.    Each  program  is  preceded  by  a  header  that  explains 

the  purpose  of  the  program,  and  describes  the  variables  passed  to  and  from  the 

program.    I  have  tried  to  make  each  header  as  complete  as  required  to  be  all  the 

documentation  necessary  to  comprehend  the  program  it  describes. 
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***************************************************************** 


* 
* 

* 
* 

** 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


Department  of  Electrical  and  Computer  Engineering 
Kansas  State  University 


* 

* 
* 

VAX  FORTRAN  source  filename:         ALIKE. FDR  * 

*************************************************************** 


ROUTINE: 


DESCRIPTION: 


DOCUMENTATION 
FILES: 


ARGUMENTS: 


X 


DEGX 


REM 


DEGR 


RETURN: 


SUBROUTINE 

ALIKE  ( I ,  J,  X,  DEGXf  REM,  DEGR) 


Refer  to  equation  (4.6)    in  the  main 
thesis.     This  program  computes  K~l_ij 
when  j=u_l,   hence  the  name  ALIKE. 


None. 


The  following  arguments  are  passed  to 
the  subroutine: 

(input)  integer 
corresponds  to  j  in  (4.6) 

( input )  intege  r 
corresponds  to  i  in  (4.6) 

(input)    real 

is  a  three  dimensional  array,  X(I,J,K) 
represents  to  Ith  coefficient  of  the 
numerator  of  the  (J,K)th  term  in  the 
partial  fraction  expansion.  Namely, 
that  term  with  the  Jth  factor  of  DEN 
to  the  Kth  power  as  denominator. 

(input)    integer 

is  an  array.     DEGX  (I,  J)    represents  the 
degree  of  the  numerator  of  the  (I,J)th 
term  in  the  partial  fraction  expansion. 
See  the  description  of  X. 

( input )    r  eal 

corresponds  to  R_l  in   (4.6)  . 

(input)    integer 

is  the  degree  of  R_l  in   (4.6)  . 

Not  used. 


ROUTINES 
CALLED: 


AUTHOR: 


PUTX 


James  F.   Stafford 


DATE  CREATED:       8Jun87     Version  1.0 


REVISIONS: 


None. 


* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
**************************************************************** 

SUBROUTINE  ALIKE  (I,  J, X,  DEGX,  REM,  DEGR) 

IMELICIT  NONE 

INTEGER  DEGR,  I,  J,  K,L, DEGX  (10,*) 

REAL*8  X(0:1,10,*),REM(0:10) 

DO  K=J,2,-1 

DEGX  (I,  K)  =DEGX  (I,  K-l) 

DO  L=0,DEGX(I,K) 

X(L,I,K)=X(L,I,K-1) 

ENDDO 
ENDDO 

CALL  HJTX  (1,1,  REM,  DEGR,  X,  DEGX) 
RETURN 
END 


ROUTINE: 


DESCRIPTION: 


SUBROUTINE 
BESSEL(F,A,N) 

This  program  computes  the  recursively 
defined  function  H_k  described  in 
chapter  V. 


DOCUMENTATION 
FILES: 


ARGUMENTS: 


***************************************************************** 

*  Department  of  Electrical  and  Computer  Engineering  * 

*  Kansas  State  University  * 

*  * 

*  VAX  FORTRAN  source  filename:         BESSEL.POR  * 
***************************************************************** 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


N 


RETURN: 


ROUTINES 
CALLED: 


AUTHOR: 


None. 


The  following  arguments  are  passed  to 
the  subroutine: 

(input)   real 

is  an  array  containing  the  coefficients 
of  the  functions  H_k.     For  a  given  j, 
F(I,0,K)   represents  the  coefficient  of 
the  cosine  term  of  H_(j-I),  with  t  to 
the  power  K.     F(I,lrK)   represents  the 
coefficient  of  the  sine  term  of  H_(j-I), 
with  t  to  the  power  K. 

( input )    r  eal 

represents  a  in  the  recursive  definition 

of  H_k  in  chapter  V. 

(input)    integer 

represents  k  in  the  recursive  definition 

of  H_k  in  chapter  V. 

Not  used. 


None. 


James  F.    Stafford 


DATE  CREATED:       9Jun37     Version  1.0 


* 

*  REVISIONS:  None. 

* 
* 

**************************************************************** 
SUBROUTINE  BESSEL(F,A,  N) 


IMPLICIT 

INTEGER 

REAL*8 

DO  J=-l,N-2 

DO  K=-2,-l 


NONE 

N,J,K 

A,  F  (-2: 0,0: 1,-1: 9) 


F(K,0,J)=F(K+1,0,J) 
F(K,1,J)=F(K+1,1,J) 


ENDDO 

ENDDO 

DO  J=0,N-2 

F(0,0,J)=F(-l,0,J)*(2*N-3)/A 
F(0,l,J)=F(-l,l,J)*(2*N-3)/A 

ENDDO 

DO  J=-l,N-3 

F(0,0fJ*-2)=F(0,0,Jl-2)-F(-2,0,J) 
F(0,l,dH-2)=F(0,l,J»-2)-F(-2,l,J) 

ENDDO 

RETURN 

END 


***************************************************************** 

*  Department  of  Electrical  and  Computer  Engineering      * 

*  Kansas  State  University  * 

*  * 

*  VAX  FORTRAN  source  filename: COMPOSE.  FOR  * 
***************************************************************** 

* 

RCUTINE:COMFOSE(X,  FACTOR,  NOM,  DEGF,  MULT) 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


DESCRIPTION: 


DOCUMENTATION 
FILES: 


ARGUMENTS: 


FACTOR 


NUM 


DEGF 


MULT 


RETURN: 


ROUTINES 
CALLED: 


AUTHOR: 


This  program  accepts  a  complex-valued 
root,X  as  input,   decides  whether  X  is 
purely  real  or  not,   and  updates  the 
factor  array,   FACTOR,    accordingly. 


None. 


(input)   complex 

Is  a  complex-valued  root  of  a  polynomial. 

( input /o  ut  put )   r  eal 

Is  an  array  containing  each  factor  of  the 

above  polynomial. 

( input /o  ut put )    int ege  r 

Is  the  number  of  factors  in  FACTOR.     NUM 

is  already  incremented  before  calling 

COMPOSE. 

( input /o  ut  put )    int  ege  r 

Is  an  array  specifying  the  degree  of  each 

corresponding  factor  in  FACTOR. 

( input/output )    integer 

Is  an  array  specifying  the  multiplicity  of 

each  factor  in  FACTOR. 

Not  used. 


None. 


James  F.    Stafford 


* 

*  DATE  CREATED:       30Jurfi8  Version  1.0 
* 

* 

*  REVISIONS:  None. 

* 

**************************************************************** 
SUBROUTINE  COMPOSE  (X,  FACTOR,  NUM,  DEGF,  MULT) 

IMPLICIT  NCNE 

INTEGER  NUM, DEGF (*)  , MULT (*) 

REAL*8  FACTOR (10, 0:2)  , SMALL 

COMPLEX*16  X 

LCGICAL  REAL 

PARAMETER  (SMALL=10E-4) 

REAL=.  FALSE. 

IF    (DREAL(X).NE.O.)    THEN 

IF    (DABS(DIMAG(X)/DREAL(X)).LT.  SMALL)    THEN 

REAL=.TRUE. 
FACTOR  ( NUM,  1)=1 
FACTOR  ( NUM,  0 )  =-DREAL  (X ) 
DEGF  (NUM)  =1 

ENDIF 

ELSE  IF    (CDABS(X).LT.  SMALL)    THEN 

REAL=.TRUE. 
FACTOR  ( NUM,  1)=1. 
FACTOR  (NUM,  0)=0. 
DEGF  (NUM)  =1 

ENDIF 

IF    (REAL.  BQ..  FALSE.)    THEN 

FACTOR  (NUM,  2 )  =1 

FACTOR  (NUM,  1 )  =-2 .  *DRE£L  (X ) 

FACTOR  ( NUM,  0 )  =OI  MG  ( X )  *DI  MAG  ( X )  +DREAL  ( X )  *DREAL  (X ) 

DEGF  (NUM)  =2 


ENDIF 

MULT(NUM)=1 
RETURN 
END 


*  Department  of  Electrical  and  Computer  Engineering      * 

*  Kansas  State  University  * 

*  * 

*  VAX  FORTRAN  source  filename:  DEFVAL.FOR  * 
***************************************************************** 

* 

ROUTINE:  COMILEX*16  FUNCTION 

DEFVAL  ( POLY,  DEG,  FACTOR,  NUMf  DEGF,  MULT,  X) 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


DESCRIPTION: 


DOCUMENTATION 
FILES: 


ARGUMENTS: 


POLY 


DEG 


FACTOR 


NUM 


DEGF 


MULT 


This  program  evaluates  a  polynomial  at 
a  given  complex  argument,  X,  all  known 
factors  are  divided  out. 


None. 


The  following  arguments  are  passed  to 
the  function: 

(input)   real 

is  an  array  containing  coefficients  of 

the  polynomial  to  be  evaluated. 

(input)    integer 

is  the  degree  of  POLY. 

(input)   real 

is  an  array  containg  the  coefficients 

of  all  known  factors  of  POLY. 

(input)    integer 

is  the  number  of  factors  in  FACTOR. 

(input)    integer 

is  an  array  specifying  the  degree  of 

each  corresponding  factor  in  FACTOR. 

(input)    integer 

is  an  array  specifying  the  multiplicity 

of  each  corresponding  factor  in  FACTOR. 

(input)  complex 

is  the  argument  at  which  the  polynomial 

is  to  be  evaluated. 


RETURN: 


ROUTINES 
CALLED: 


AUTHOR: 


DATE  CREATED: 


REVISIONS: 


Not  used. 

None. 

James  F.    Stafford 

30Jun38  Version  1.0 

None. 


* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
***************************************************************** 

COMILEX*16  FUNCTION  DEFVAL  (POLY,  DEC,  FACTOR,  NUM,  DEGF,  MULT,  X) 


NONE 

DEC,  I,  NUM,  DEGF  (* )  , MULT (* ] 

POLY(0:*) , FACTOR (10, 0:2) 

X,EVAL 


IMPLICIT 

INTEGER 

REAL*8 

COMPLEX*16 

DEFVAL=POLY(DEG) 

DO   I=DEG-1,0,-1 

DEEVAL=DEFVAL*X+POLY  (I ) 
ENDDO 
DO  1=1, NUM 

IF    (DEFVAL. NE.0)    THEN 

DEFVAL=DEFVAL/EVAL(  FACTOR,  I,  DEGF  (I)  ,MULT(I)  ,X) 

ENDIF 
ENDDO 
RETURN 
END 
COMPLEX*16  FUNCTION  EVAL  (FACTOR,  I,  DEGF,  MULT,  X) 


IMPLICIT  NONE 

INTEGER  J,  I,  DEGF,  MULT 

REAL*8  FACTOR(10,0:2) 

G0MHLEX*16  Xf  VALUE 

EVAL=1 
VALUE=FACTOR  ( I ,  DEGF) 

DO  J=DEGF-1,0,-1 

VAL  UE=VAL  UE*X+FACTOR  ( I ,  J) 
ENDDO 
DO  J=1,MULT 

EVAL=EVAL*VALUE 
ENDDO 
RETURN 
END 


***************************************************************** 


*  Department  of  Electrical  and  Computer  Engineering      * 

*  Kansas  State  University  * 

*  * 

*  VAX  FORTRAN  source  f  ilename:DER]V.POR  * 
***************************************************************** 

* 

*  ROUTINE:  SUBROUTINE 

DERIV(  POLY,  DEC) 


* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


DESCRIPTION: 


DOCUMENTATION 
FILES: 


ARGUMENTS: 


POLY 


DEC 


RETURN: 


ROUTINES 
CALLED: 


AUTHOR: 


DATE  CREATED: 


REVISIONS: 


This  program  computes  the  derivative  of 
a  given  polynomial. 


None. 


The  following  arguments  are  passed  to 
the  routine. 

( input /o  ut  put )   r  eal 

is  an  array  containing  the  coefficients 
of  the  polynomial  to  be  differentiated. 
On  return,   this  array  contains  the 
coefficients  of  the  derivative. 

( input /out put )    intege r 
is  the  degree  of  the  polynomial  on 
input  and  the  degree  of  the 
derivative  on  output. 


Not  used. 


None. 


James  F.    Stafford 


30Jun88  Version  1.0 


None. 


**************************************************************** 
SUBROUTINE  DEPJV(POLY,DEG) 


IMPLICIT 

NONE 

INTEGER 

I,  J,EEG 

REAL*8 

POLY(0:*) 

DO  J=0,DEG-1 

POLY  (J)  = 

(JH)*POLY(JH) 

ENDDO 

POLY(DEG)=0. 

DEG=DEG-1 

RETURN 

END 

***************************************************************** 
*  Department  of  Electrical  and  Computer  Engineering  * 


Kansas  State  University 


* 

*  * 

*  VAX  FORTRAN  source  filename:         DIFFERENT.  FOR  * 
***************************************************************** 

* 

ROUTINE:  SUBROUTINE 

DIFFERENT  (I,  B,  DEGB,  DEN,  DEGD,  MULTS,  X,  DEGX) 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


DESCRIPTION: 


DOCUMENTATION 
FILES: 


ARGUMENTS: 


B 


DEGB 


DEN 


DEGD 


MULTS 


Refer  to  equation  (4.6)    in  the  main 
thesis.     This  program  computes  K~l_ij 
when  j<u_l,   hence  the  name  DIFFERENT. 


None. 


The  following  arguments  are  passed  to  the 
subroutine: 

(input)    integer 
corresponds  to  j   in  (4.6) 

(input)    real 

is  an  array  containing  the  coefficients  of 
the  polynomial  f_l  in  the  notation  of 
chapter  IV. 

( input )    intege  r 
is  the  degree  of  B 

( input )    r  eal 

is  a  two-dimensional  array.     DEN  (I,  J) 
represents  the  coefficient  of  the  Ith 
power  of  x  in  the  Jth  factor  of  the 
denominator  polynomial. 

(input)    integer 

is  an  array.     DEGD(I)    represents  the 
degree  of  the  Ith  factor  in  the 
denominator  polynomial. 

(input)    integer 

is  an  array.  MULTS(I)  represents  the 
multiplicity  of  the  Ith  factor  in  the 
denominator  polynomial. 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

** 


DEGX 


RETURN: 


ROUTINES 
CALLED: 


AUTHOR: 


(input)   real 

is  a  three-dimensional  array.     X(I,J, K) 
represents  the  Jth  coefficient  of  the 
numerator  of  the  (JfK)th  term  in  the 
partial  fraction  expansion.     Namely, 
that  term  with  the  Jth  factor  of  DEN 
to  the  Kth  power  as  denominator. 

(input)    integer 

is  a  two-dimensional  array.     DEG(I,J) 
represents  the  degree  of  the  numerator 
of  the   (I,J)th  term  in  the  partial 
fraction  expansion.     See  the  description 
of  X. 

Not  used. 


EUCLID,  HJTX,  GETD,  GETX,  POLADD 


James  F.  Stafford 


DATE  CREATED:       8Jun37     Version  1.0 


REVISIONS: 


None. 


************************************************************** 
SUBROUTINE  DIFFERENT  (I,  BfDEGB,  DEN,  DEGD,MULTS,  X,  DEGX) 


IMFLICIT 
INTEGER 

REAL*8 


NONE 

DEGD  ( * )  ,  MULTS  ( * )  ,  I,  J,  K,  L,  DECS, 
DEGT,DEGX(10,*)  ,DEGA,DEGB,DEGF 

DEN(0:2,*)  ,X(0:1,10,*)  ,A(0:2)  ,B(0:2)  , 
S(0:1)  ,T(0:1)  ,F(0:1) 


DO  J= 1-1,1,-1 

PRINT  *,  ,J=I,J 

CALL  GETD  (J,  A,  DEGA,  DEN,  DEGD) 

CALL  GETX(J,MULTS(J)  ,F,DEGF,X,DEGX) 


DO  K=MULTS( J) -1,1,-1 
PRINT  VK=',K 

CALL  EUCLID  (A,  DEGA,  B,  DEGB,  F,  DEGF,  S,  DEGS,  T,  DEGT) 
CALL  PUTX(J,K+1,T,DEGT,X,DEGX) 
CALL  GETX(J,K,F,DEGF,X,DEGX) 
CALL  POLADD(F,  DEGF,  S,  DEGS) 
ENDDO 

CALL  EUCLID  (A,  DEGA,  B,  DEGB,  F,  DEGF,  S,  DEGS,  T,  DEGT) 
CALL  PUTX(J,1,T,DEGT,X,DEGX) 
CALL  GETX(I,1,F,DEGF,X,DEGX) 
CALL  POLADD (F,  DEGF,  S,  DEGS) 
CALL  PUTX(I,1,F,DEGF,X,DEGX) 
ENDDO 
RETURN 
END 


ROUTINE: 


SUBROUTINE 

EUCLID (POL1 ,DEG1 , POL2 ,DEG2 , GCD, DEGG) 


DESCRIPTION:         This  program  computes  the  greatest 

common  divisor  of  two  given  polynomials. 


DOCUMENTATION 

FILES:  None. 


***************************************************************** 

*  Department  of  Electrical  and  Computer  Engineering  * 

*  Kansas  State  University  * 

*  * 

*  VAX  FORTRAN  source  filename: EUCLID.  FOR  * 
***************************************************************** 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


ARGUMENTS: 
P0L1 

DEG1 
POL2 

DEG2 
GCD 


DEGG 


RETURN: 


ROUTINES 
CALLED: 


(input)   real 

is  an  array  representing  one  input 

polynomial. 

(input)    integer 

is  the  degree  of  POLL 

(input)   real 

is  an  array  representing  the  other 

input  polynomial. 

(input)    integer 

is  the  degree  of  POL2. 

(output)   real 

is  the  gcd  of  the  two  input  polynomials 

(output)  integer 

is  the  degree  of  GCD. 


Not  used. 


AUTHOR: 


POLDIV 


James  F.  Stafford 


* 
* 

*  DATE  CREATED:   30JunS8  Version  1.0 

* 
* 

*  REVISIONS:  None. 
* 

* 
**************************************************************** 

SUBROUTINE  EUCLID (POLlfDEGlfPOL2/DEG2,GCD,DEGG) 

IMPLICIT  NCNE 

INTEGER  If  J,DEG1,DEG2,DEGG,DEGA,DEGB,DEGQ 

REAL*8  POL1(0:10)  ,POL2(0:10)  ,GCD(0:10)  , 

+  A(0:10)  ,B(0:10)  ,Q(0:10)  ,ZERO 

LOGICAL  EASY,  HARD 

PARAMETER  (ZERO=l  .OE-5) 

E  AS  Y=.  FALSE. 
HARD=.TRUE. 

CALL  GET(GCD,DEGG,POLl,DEGl) 
CALL  GET(B,DEGB,POL2,DEG2) 

DO  WHILE    (.NOT.  (DEGB.  EQ.O  .AND.DABS(B(0)  ).LT.ZERO)) 

CALL  GET(A,DEGA,GCD,DEGG) 

CALL  GET (GCD, DEGG,B, DEGB) 

CALL  POLDIV  (A,  DEGA,  GCD,  DEGG,Q,DEGQ,B,  DEGB,  EASY,  HARD) 

DO   I=0,DEGG 

GCD  (I )  =GCD  (I ) /GCD  (DEGG) 

ENDDO 
ENDDO 
RETURN 
END 
SUBROUTINE  GET  (A,  DEGA,  B,  DEGB) 


IMHilCIT  NCNE 

INTEGER  I ,  J,  DEGA,  DEGB 

REAL*8  A(0:*),B(0:*) 

DEGA=DEGB 

DO  1=0,  DEGA 

A(I)=B(I) 
ENDDO 
RETURN 
END 


ROUTINE: 


DESCRIPTION: 


SUBROUTINE 

EUCLID  (A,  DEGAf  B,  DEGB,  F,  DEGF,  S,  DEGS, 

TfDEGT) 

This  program  computes  A  and  A"*  in 
equation   (4.1)    given  Kf    s  and  r  via 
the  methods  discussed  in  chapter  IV. 


DOCUMENTATION 
FILES: 


ARGUMENTS: 


None. 


***************************************************************** 

*  Department  of  Electrical  and  Computer  Engineering  * 

*  Kansas  State  University  * 

*  * 

*  VAX  FORTRAN  source  filename:  EUCLIDEAN.  FOR  * 
***************************************************************** 
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* 

* 


DEGA 


B 


DEGB 


DEGF 


DEGS 


The  following  arguments  are  passed  to 
the  subroutine: 

( input )    r  eal 

is  an  array  containing  the  coefficients 

of  s  in   (4.2) . 

(input)    integer 

is  the  degree  of  s  in   (4.2) . 

(input)   real 

is  an  array  containing  the  coefficients 

of  r  in   (4.2) . 

(input)    integer 

is  the  degree  of  r  in   (4.2) . 

(input)   real 

is  an  array  containing  the  coefficients 

of  K  in   (4.2) . 

(input)    integer 

is  the  degree  of  K  in   (4.2). 

(output)   real 

is  an  array  containing  the  coefficients 

of  A  in   (4.2) . 

(output)    integer 

is  the  degree  of  A  in   (4.2)  . 


* 
* 

* 
* 
* 
* 
* 
* 
* 
* 
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* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
**************************************************************** 

SUBROUTINE  EUCLID  (A,  DEGA,  B,  DEGB,  F,  DEGF,  S,  DEGS, 

+  T,  DEGT) 


DEGT 


RETURN: 


ROUTINES 
CALLED: 


AUTHOR: 


DATE  CREATED: 


REVISIONS: 


(output)   real 

is  an  array  containing  the  coefficients 

of  A"*  in   (4.2)  . 

(output)    integer 

is  the  degree  of  A~*  in  (4.2) . 

Not  used. 


POLDIV,K)LMULT, 
GET(contained  in  TRANSFER) 


James  F.    Stafford 
9Jun87     Version  1.0 
None. 


IMILICIT 
INTEGER 

REAL*8 

LOGICAL 

ADD=0 

DIV=0 

MULT=0 

EASYDIV=.TRUE. 

HARD=.  FALSE. 

DEGR=1 

IF    (DEGA.  LE.  DEGB)      THEN 

SWITCH=.TRUE. 


NONE 

I ,  DEGA,  DEGB,  DEGF,  DEGT,  DEGS,  DEGQ,  DEGR, 
DEG1 ,  DEG2 ,  ADD,  MULT,  DIV 

A(0:2)  ,B(0:2)  ,F(0:1)  ,S(0:1)  ,T(0:1)  , 
QUO(0:1)  ,REM(0:2) ,BUFF1(0:3)  ,BUFF2(0:3) 

EASYDIV,  HARD,  EASYMULT,  SWITCH 


CALL  GET(BUFF2,DEG2,A,DEGA) 
CALL  GET(BUFF1,DEG1,B,DEGB) 

ELSE 

SWITCH=.  FALSE. 

CALL  GET(BUFF1,DEG1,A,DEGA) 
CALL  GET(BUFF2,DEG2,B,DEGB) 

ENDIF 

1=0 

DO  WHILE    (DEGR.GT.0) 

CALL  POLDIV  ( BUFFI ,  DEGl ,  BUFF2 ,  DEG2 ,  QUO,  DEGQ ,  REM,  DEGR, 
EASYDIV,  HARD) 

CALL  GET  (BUFFI,  DEGl  ,BUFF2,DEG2) 
CALL  GET(BUFF2,DEG2f REM, DEGR) 

E  AS  YDIV=.  FALSE. 

HARD=.TRUE. 

1=1+1 

ENDDO 

IF    (I.LT.2)   THEN 

EASYMULT=.TRUE. 
DEGQ=0 

ELSE 

EASYMULT=.  FALSE. 

BUFF2(0)=-BUFF2(0) 

ADD=ADD«-1 

ENDIF 

CALL  POLMULT( QUO,  DEGQ ,F,DEGF,  BUFFI, DEGl  rEASYMULT) 

DO   1=0,  DEGl 

BUFFI ( I ) =BUFF1 ( I ) /BUFF2 ( 0) 
DIV=DIV+1 


ENDDO 


HARD=.  FALSE. 

IF    (SWITCH)   THEN 

CALL  POLDIV(BUFFl,DEGl,A,DEGA,Cro,DEC^,T,DEGT, 
EASYDIV,HARD) 

DEGS=DEGB-1 
S(DEGS)=T(DEGT) 
DO  I=DEGS-1,0,-1 

S  ( 0 )  =T  (DEGT)  *  ( B  (DEGB-1 )  -A  (DEGA-1 )  ) 

MULT=MULT+1 

ADD=ADDfl 

IF    (DEGT.GT.O)    THEN 

S(0)=S(0)+T(0) 
ADD=ADDfl 

ENDIF 

ENDDO 

ELSE 

CALL  POLDIV( BUFFI ,DEGl,Bf  DEGB,  QUO, DEGQ,S,DEGS, 
EASYDIV,HARD) 

DEGT=1 

S(0)=>-S(0) 

T(1)=S(0) 

T(0)=(A(1)-B(0))*S(0)+F(1) 

ADD=ADEH-3 

MULT=MULT+1 

ENDIF 

PRINT  *, ADD,  'additions1 
PRINT  *,MQLT,  'multiplies' 
PRINT  *,DIV,  'divisions' 

RETURN 

END 
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Department  of  Electrical  and  Computer  Engineering  * 

Kansas  State  University  * 

* 

VAX  FORTRAN  source  filename:    EXPAND. FOR  * 

*************************************************************** 


ROUTINE: 


DESCRIPTION: 


DOCUMENTATION 
FILES: 


ARGUMENTS: 


NUM 


DEGN 


DEN 


DEGD 


MULTS 


SUBROUTINE 

EXPAND ( (NUM,  DEGN,  DEN,  DEGD,  MULTS,  X,DEGX, 

NO_  FACTS) 


This  program  performs  a  partial  fraction 
expansion  on  a  rational  function  using 
Chin's  and  Steiglitz's  algorithm. 


None. 


The  following  arguments  are  passed  to 
the  subroutine: 

(input)   real 

is  an  array  containing  the  coefficients 
of  the  numerator  polynomial  of  the 
rational  function  to  be  expanded. 

(input)  integer 

is  the  degree  of  the  numerator 

polynomial. 

(input)   real 

is  a  two-dimensional  array.     DEN (I, J) 
represents  the  coefficient  of  the  Ith 
power  of  x  in  the  Jth  factor  of  the 
denominator  polynomial. 

(input)    integer 

is  an  array.     DEGD(I)   represents  the 
degree  of  the  Ith  factor  in  the 
denominator  polynomial. 

(input)    integer 

is  an  array.  MULTS  (I)  represents  the 
multiplicity  of  the  Ith  factor  in  the 
denominator  polynomial. 


NO_ FACTS         (input)    integer 
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DEGX 


RETURN: 


ROUTINES 
CALLED: 


AUTHOR: 


is  the  number  of  factors  in  the 
denominator  polynomial. 

(output)   real 

is  a  three-dimensional  array.     X(If  JfK) 
represents  the  Ith  coefficient  of  the 
numerator  of  the  (JfK)th  term  in  the 
partial  fraction  expansion.     Namely, 
that  term  with  the  Jth  factor  of  DEN 
to  the  Kth  power  as  denominator. 

(output)    integer 

is  an  array.     DEGX(I,  J)    represents  the 
degree  of  the  numerator  of  the  (I,  J)th 
term  in  the  partial  fraction  expansion. 
See  the  description  of  X. 


Not  used. 


GETD,  POLDIV ,  ALIKE ,  DI FFERENT 


James  F.    Stafford 


DATE  CREATED:       6JunS7     Version  1.0 


REVISIONS: 


None. 


************************************************************* 

SUBROUTINE  EXPAND  (NUMf  DEGN,  DEN,  DEGD,  MULTS,  X,  DEGX, 

+  NO_  FACTS) 


+ 
+ 


IMPLICIT 
INTEGER 

REAL*8 

LOGICAL 
EASY=.  FALSE. 


NONE 

DEGN,  NO_  FACTS,  DEGD(  10)  ,DEGB, 
MULTS  ( * )  , DEGR,  I ,  DEGQ ,  J,  K,  L, 
DEGX(10,*) 

NUM(0:*) ,DEN(0:2,*),X(0:1,10,*) , 
QUO(0:10)  ,REM(0:10)  ,B(0:2) 

EASY,  HARD 


HARD=.  FALSE. 
DO  1=1  ,NQ_  FACTS 

PRINT  *,'I=', I 

CALL  GETD(I,B,DEGB,DEN,DEGD) 

DO  J=1,MJLTS(I) 

CALL  POLDW(NUM,DEGN,B,DEGB,CTO,DEGQ,REM,DEGR, 
EASY,  HARD) 

DEGN=DEGQ 

DO  K=0,DEGN 

NUM(K)=OUO(K) 

ENDDO 

PRINT  *,DEGRf,YESl 

DO  K=0,DEGR 

PRINT  *,REM(K) 

ENDDO 

CALL  ALIKE(I,J,X,DEGX,REM,DEGR) 

CALL  DIFFERENT(I,B,DE£B,DEN,DE£D,MULTS,X,DEGX) 

ENDDO 

ENDDO 

RETURN 

END 


***************************************************************** 


Department  of  Electrical  and  Computer  Engineering      * 
Kansas  State  University  * 

* 

VAX  FORTRAN  source  f  ilename:FACTORER.FOR  * 

***************************************************************** 

* 

*  ROUTINE: 

* 

* 

* 


SUBROUTINE 

FACTORER  (POLY,  DEGP,  FACTOR,  NUM,  DEGF,  MULT) 


* 

* 
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* 
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* 

* 
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DESCRIPTION: 


DOCUMENTATION 
FILES: 


ARGUMENTS: 
POLY 

DEGP 
FACTOR 

NUM 
DEGF 

MULT 

RETURN: 


ROUTINES 
CALLED: 


This  program  factors  a  given  input 
polynomial  into  irreducible  elements  of 
R[x]. 


None. 


( input )   r  eal 

is  an  array  containing  the  coefficients 

the  polynomial  to  be  factored. 

(input)    integer 

is  the  degree  of  POLY. 

(output)   real 

is  an  array  containinng  coefficients 

each  factor  of  POLY. 

(output)    integer 

is  the  number  of  factors  in  FACTOR. 

(output)    integer 

is  an  array  specifying  the  degree  of 

the  corresponding  factor  in  FACTOR. 

(output)  integer 

is  an  array  specifying  the  multiplicity 

of  the  corresponding  factor  in  FACTOR. 

Not  used. 


DERIV,    EUCLID,    FIND_EM 


* 

*  AUTHOR:        James  F.  Stafford 

* 
* 

*  DATE  ORFATED:        30Jun88   Version  1.0 

* 
* 

*  REVISIONS:  None. 
* 

* 
************************************************************* 

SUBROUTINE  FACTORER  (POLY,  DEGP,  FACTOR,  NUM,  DEGF,  MULT) 

IMPLICIT  NONE 

INTEGER  I,J,K,DEGP,DEGF(*),MUXT(*),NUM,DEGGCD(0:10)  , 

+  DEGD, DEGG 

REAL*8  FOLY(0:10)  ,FACTOR(10,0:2)  ,GCD(0:10,0:10)  ,D(0:10) 

+  G(0:10) 

DO  1=0,  DEGP 

D(I)=POLY(I) 
GCD(0fI)=POLY(I) 

ENDDO 

DEGD=DEGP 

DEGGCD(0)=DEGP 

K=0 

DO  WHILE    (DEGGCD(K).GT.O) 

K=K+1 

CALL  DERIV(D,DEGD) 

CALL  EUCLID  (POLY,  DEGP,  D,  DEGD,  G,  DEGG) 

DO  J=0,DEGG 

GCD(K,J)=G(J) 

ENDDO 

DEGGCD(K)=DEGG 
ENDDO 
DO   I=K-1,0,-1 


DO  J=0,DEGGCD(I) 
G(J)=GCD(I,J) 

ENDDO 

DEGG=DEGGCD(I) 

CALL  FIND. EM ( G,  DEGG ,  FACTOR,  NUM,  DEG F,  MULT) 
ENDDO 
RETURN 
END 


***************************************************************** 

*  Department  of  Electrical  and  Computer  Engineering  * 

*  Kansas  State  University  * 

*  * 

*  VAX  FORTRAN  source  filename: FIND. EM.  FOR  * 
***************************************************************** 

* 

SUBROUTINE 

FIND_  EM  (P,  DEGPf  FACTOR,  NUMf  DEGF,  MULT) 


ROUTINE: 


DESCRIPTION: 


DOCUMENTATION 
FILES: 


ARGUMENTS: 


RETURN: 


ROUTINES 
CALLED: 


AUTHOR: 


DATE  CREATED: 


REVISIONS: 


A  polynomial  and  a  list  of  already  known 
roots  is  passed  to  the  subroutine. 


None. 
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**************************************************************** 

SUBROUTINE  FIND_EM(P,  DEGPf  FACTOR,  NUM,  DEGF ,  MULT) 


Not   used. 


MULLER 


James  F.    Stafford 
5Sep86     Version  1.0 
None. 


IMPLICIT 

INTEGER 

REAL*8 

SUM=0 

DO  J=lfNUM 


NONE 

J,DEGP,NUM,DEGF(*)  ,MULT(*)  ,SUM 

P(0:*)  ,FACTOR(10f0:2) 


MULT(J)=MULT(J)+1 
SUM=SUW-DEGF(J)  *MULT(  J) 

ENDDO 

DO  WHILE    (DEGP-SUM.GT.O) 

CALL  MULLER  (P,  DEGP,  FACTOR,  NUM,  DEGF,  MULT) 
SUM=  SUM- DEGF  (NUM) 

ENDDO 

RETURN 

END 


***************************************************************** 


*  Department  of  Electrical  and  Computer  Engineering  * 

*  Kansas  State  University  * 

*  * 

*  VAX  FORTRAN  source  filename:         GAMMA.  FOR  * 
***************************************************************** 

* 

FUNCTION 
GAMMA  (K) 
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* 

******* 


ROUTINE: 


DESCRIPTION: 


DOCUMENTATION 
FILES: 


ARGUMENTS: 


K 


RETURN: 


ROUTINES 
CALLED: 


AUTHOR: 


This  program  computes  the  integer-valued 
function  Gamma  defined  in  the  glossary. 


None. 


The  following  argument  is  passed  to  the 
function: 

(input)    integer 

is  any  non-negative  integer. 

The  function  returns  an  integer  according 
to  the  definition  in  the  glossary. 


None. 


James  F.  Stafford 


DATE  CREATED:        5Sep36     Version  1.0 


REVISIONS: 


None. 


********************************************************* 


INTEGER  FUNCTION 
IMPLICIT  NONE 

INTEGER  J,K 

GAMMA=1 


GAMMA(K) 


DO  J=K-1,2,-1 

GAMMA=GAMMA*J 
ENDDO 
REIUEN 
EM) 


***************************************************************** 

*  Department  of  Electrical  and  Computer  Engineering      * 

*  Kansas  State  University  * 

*  * 

*  VAX  FORTRAN  source  filename:         INIT.  FOR  * 
***************************************************************** 

* 

ROUTINE:  SUBROUTINE 

INIT(F) 


* 
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**** 


DESCRIPTION: 


DOCUMENTATION 
FILES: 


ARGUMENTS: 


RETURN: 


ROUTINES 
CALLED: 


AUTHOR: 


This  program  initializes  the  recursively 
defined  function  H_k  described  in 
chapter  V. 


None. 


The  following  argument  in  passed  to 
the  subroutine: 

(output)   real 

is  an  array  containing  the  coefficients 
of  the  functions  H_k.     For  a  given  j, 
F(I,0,K)    represents  the  coefficient  of 
the  cosine  term  of  H_(j-I),  with  t  to  the 
power  K.     F(I,1,K)   represents  the 
coefficient  of  the  sine  term  of  H_(j-I) 
with  t  to  the  power  K. 

Not  used. 


None. 


James  F.    Stafford 


DATE  CREATED:        9Jun37     Version  1.0 


REVISIONS: 


None. 


************************************************************ 
SUBROUTINE  INIT(F) 


IMILICIT  NONE 

INTEGER  I,J,K 

REAL*8  F(-2:0, 0:1,-1:9) 

DO  I=-2,0 

DO  J=0,1 

DO  K=-l,9 

F(I,J,K)=0. 

ENDDO 

ENDDO 

ENDDO 

F(-1,0,-1)=1. 
F(0,1,0)=1. 

RETURN 

END 


DESCRIPTION: 


***************************************************************** 

*  Department  of  Electrical  and  Computer  Engineering      * 

*  Kansas  State  University  * 

*  * 

*  VAX  FORTRAN  source  filename:  INPUT. RAT.  FOR  * 
***************************************************************** 

* 

*  ROUTINE:  PROGRAM 
* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

**************************************************************** 

NONE 


DOCUMENTATION 
FILES: 


ARGUMENTS: 


RETURN: 


ROUTINES 
CALLED: 


AUTHOR: 


DATE  CREATED: 


REVISIONS: 


This  program  allows  one  to  establish  a 
data  file  compatible  with  the  inverse 
transform  package  programs  containing  the 
necessary  data  to  describe  a  rational 
function. 


None. 
Not  used. 
Not  used. 

None. 

James  F.  Stafford 

27May87  Version  1.0 

None. 


IMPLICIT 
INTEGER 
REAL*8 
COMPLEX*16 
CHARACTER* 15 


I ,  N_  DEG ,  D_  DEG ,  NO_  ROOTS ,  MLTPLCTS  ( 1 0) 
NUM(0:15)  ,DEN(0:15) 
ROOTS  (10) 
FILENAME, YESNO 


PRINT  *,'Are  numerator  roots  known?   (Y/N)  • 

READ   (*,200)    YESNO 

IF    (YESNO.  EQ.'Y')    THEN 

CALL  INPUT. FACT  ( NO. ROOTS,  ROOTS,  MLTFLCTS ) 

CALL  RECONSTRUCT  (NO_ ROOTS,  ROOTS,  MLTFLCTS, 
+  NUM,N_DEG) 

ELSE 

CALL  rNFUT_NONFACT(NUM,N.DEG) 
ENDIF 

PRINT  *,'Are  denominator  roots  known?   (Y/N)' 
READ    (*,200)    YESNO 
IF    (YESNO.  EQ.'Y')    THEN 

CALL   INPUT. FACT  (NO_ ROOTS,  ROOTS,  MLTPLCTS) 

CALL  RECONSTRUCT  (NO.  ROOTS,  ROOTS,  MLTPLCTS, 
+  DEN,D_DEG) 

ELSE 

CALL  INPUT_NCNFACr(DEN,D_DEG) 

ENDIF 

PRINT  *,' Enter  filename.' 

READ   (*,200)    FILENAME 

200  FORMAT    (A15) 

OPEN   (UNIT=1,    FILE=FILENAME,    STATUS='NEW') 

WRITE    (1,*)    N.DEG 

DO  1=0,  N_  DEC 

WRITE    (1,*)    NUM(I) 

ENDDO 


WRITE    (1,*)   D_DEG 

DO  I=0,D_DEG 

WRITE    (1,*)   DEN(I) 

ENDDO 

CLOSE   (UNIT=1,   STATUS^KEEP' ) 

END 

SUBROUTINE  INFUT_FACT  (NO_ ROOTS,  ROOTS,  MLTPLCTS) 

IMELICIT  NONE 

INTEGER  NO_  ROOTS,  MLTPLCTS  (*),  I 

COMPLEX*16  ROOTS  (*) 

PRINT  *, '  Input  nimber  of   roots' 

READ   (*,*)    NO_ ROOTS 

100  FORMAT    (F8.5,F8.5) 

DO  1=1, NO_  ROOTS 

PRINT  *,' Input  root  nimber  ',1 

READ    (*,100)   ROOTS(I) 

PRINT  *,' Input  corresponding  multiplicity' 

READ    (*,*)   MLTPLCTS(I) 

PRINT  *,ROOTS(I),MLTPLCTS(I) 

ENDDO 

RETURN 

END 

SUBROUTINE  RECONSTRUCT  (NO. ROOTS,  ROOTS,  MLTPLCTS, 

+  RESULT,  ORDER) 

IMPLICIT  NONE 

INTEGER  NQ_ ROOTS,  MLTPLCTS (*), I,  J,  ORDER,  MULT 

COMPLEX*16  ROOTS (*)  , BUFFER (0:50) 


REAL*8  RESULT(0:*) 

BUFFER  ( 0)  =DCMPLX  ( 1 . 0 , 0 . 0 ) 

DO  1=1,50 

BUFFER  ( I )  =DCMHJX  ( 0 . ,  0  . ) 

ENDDO 

ORBER=0 

DO   1=1  ,NO_  ROOTS 

ORDER=ORDER+  MLTPLCTS  ( I ) 
MULT=MLTPLCTS(I) 

DO  WHILE    (MULT.NE.O) 

DO  J=GRDER,1,-1 

BUFFER  ( J)  =BUFFER  ( J-l )  -BUFFER  ( J)  *ROOTS  ( I ) 

ENDDO 

BUFFER ( 0)  =-BUFFER  (0)  *ROOTS  ( I ) 
MULT=MULT-1 

ENDDO 
ENDDO 
DO  1=0,  ORDER 

RESULT  ( I )  =REAL  ( BUFFER  ( I ) ) 

ENDDO 

RETURN 

END 

SUBROUTINE  INHJT_  NONFACT  (POLY,  DEG) 

IMPLICIT  NONE 

INTEGER  DEG,  I 

REAL*8  POLY(0:*) 


PRINT  *, '  Input  degree' 

READ  (*,*)  DBG 

DO  1=0, DEG 

PRINT  *,  'Input  coeff.   of  power  '  ,1 
READ   (*,*)    POLY  (I) 

ENDDO 

RETURN 

END 


ROUTINE: 


DESCRIPTION: 


***************************************************************** 

*  Department  of  Electrical  and  Computer  Engineering  * 

*  Kansas  State  University  * 

*  * 

*  VAX  FORTRAN  source  filename:         INV2.F0R  * 
***************************************************************** 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


DOCUMENTATION 
FILES: 


ARGUMENTS: 


SUBROUTINE 

INV2  ( I ,  J,  RESP,  OMEGA,  A,  B) 

This  program  computes  the  inverse 
Laplace  transform  of  a  rational  function 
such  that  equation  (5.2)   applies.     Note 
that  since  the  algorithm  described  in 
chapter  V  is  recursive,   that  for  a 
given  kf   all  of  the  previous  transforms 
for  k>j>=l  must  be  already  computed. 
The  function  H_k  is  computed  each  time 
the  subroutine  is  called,   using  H_(k-1) 
and  H_  (k-2)   which  are  held  in  an  array 
intrinsic  to  this  routine,   namely,   F. 
On  K=l,   F  is  initialized  to  hold  the 
coefficients  of  H_0  and  H_-l.  The 
inverse  transform  coefficients  are 
accumulated  in  an  array  called  RESP. 


None. 


RESP 


OMEGA 


The  following  arguments  are  passed  to 
the  subroutine: 

(input)    integer 

is  an  index  variable  specifying  which 
factor  of  the  denominator  polynomial 
of  the  original  rational  function  is 
currently  being  inverse  transformed. 

(input)    integer 
corresponds  to  k  in  (5.2) 

( input /o  ut  put )   r  eal  \ 
is  an  array  to  accumulate  the  computed 
time-domain  response.     RESP(j,0,i) 
represents  alpha_ji  in  equation  (5.5) 
and  RESP(j,l,i)    respresents  beta_ji 
in  equation  (5.5)  . 

(input)   real 


B 


RETURN: 


ROUTINES 
CALLED: 


AUTHOR: 


DATE  CREATED: 


REVISIONS: 


corresponds  to  a  in  (5.2) . 

(input)   real 

corresponds  to  A  in   (5.2) . 

(input)  real 

corresponds  to  B  in  (5.2) . 

Not  used. 


INIT,  BESSELf  GAMMA 
James  F.   Stafford 
5Sep36     Version  1.0 
None. 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

**************************************************************** 

SUBROUTINE  INV2  (I,  J,RESP,  OMEGA,  A,B) 

IMPLICIT  NONE 

INTEGER  I,J,K,L,M,GAMMA 

REAL*8  A, B, OMEGA, F (-2: 0,0: 1,-1: 9)  ,RESP(10,0:1,0:9)  , 

+  ADJ 

IF    (J.EQ.l)      THEN 

CALL   INIT(F) 
ELSE 

CALL  BESSEL(F,OMEGA,J) 

ENDIF 

ADJ=GAMMA(  J)  *  (2*OMEGA)**  (J-l) 

DO  K=0,J-1 

RESP(I,0,K)=RESP(I,0,K)+(F(-1,0,K-1)*A+F(0,0,K)*B) 
+  /ADJ 


RESP(I,l,K)=RESP(I,l,K)+(P(-lrl,K-l)*AfF(0,l/K)*B) 

+  /AECT 

ENDDO 

RETURN 

END 


***************************************************************** 


*  Department  of  Electrical  and  Computer  Engineering      * 

*  Kansas  State  University  * 

*  * 

*  VAX  FORTRAN  source  filename:         INVERT. FOR  * 
***************************************************************** 

* 

PROGRAM 

INVERT 


ROUTINE: 


DESCRIPTION: 


DOCUMENTATION 
FILES: 


ARGUMENTS: 


This  program  computes  the  inverse 
Laplace  transform  of  a  rational 
function  that  has  been  expanded  into 
partial  fractions.     The  user  is  prompted 
for  a  filename  under  which  the 
output  of  the  partial  fraction  expander 
program  has  been  stored.     Again,   the 
user  is  prompted  for  another  filename 
under  which  to  store  the  parameters  of 
the  inverse  transform  function. 


None. 


RETURN: 


ROUTINES 
CALLED: 


AUTHOR: 


None. 


Not  used. 


INV2 


James  F.  Stafford 


DATE  CREATED:        10JunS7  Version  1.0 


REVISIONS: 


None. 


**************************************************************** 
PROGRAM  INVERT 


IMPLICIT 
INTEGER 


NONE 


I,  J,  K,NO_ TERMS,  ORDER,  MULT(IO)  , GAMMA 


REAL*8  TAU(IO)  ,OMEGA(10)  ,A,B,F(-2: 0,0:1,-1:10)  , 

+  RESP(10,0:1,0:9) 

CHARACTER* 15         FILENAME 

PRINT  *,  *  Enter  filename. ' 

READ   (*,200)    FILENAME 

200  FORMAT    (A15) 

OPEN   (UNIT=1,    FILE=FILENAME,    STATU S= '  OLD1 ) 

READ   (1,*)    NO_ TERMS 

DO  1=1,  NO_  TERMS 

READ   (1,*)   ORDER 
PRINT  *,  ORDER 
READ    (1,*)    MULT(I) 
PRINT  *,MULT(I) 

IF    (ORDER.  EQ.l)    THEN 

READ    (1,*)   TAU(I) 
PRINT  *,TAU(I) 

ELSE 

READ   (1,*)   TAU(I) 
PRINT  *,TAU(I) 
READ   (1,*)   OMEGA(I) 
PRINT  *,OMEGA(I) 

ENDIF 

DO  J=1,MULT(I) 

IF    (ORDER.  EQ.l)    THEN 

READ   (1,*)    A 

PRINT  *,A 

RESP  ( 1 , 0 ,  J-l )  =A/GAMMA  ( J) 


ELSE 


READ  (1,*)  A 
PRINT  *,'A  =',A 
READ  (1,*)  B 


PRINT  *,  'B  =',B 

CALL  INV2(I,J,RESP,0MEGA(I),A,B) 

ENDIF 

ENDDO 

ENDDO 

CLOSE    (UNrT=l,    STATU S=' KEEP' ) 

PRINT  *,'  Enter  filename.1 

READ    (*,200)    FILENAME 

OPEN   (UNIT=1,   FILE=FILENAME,    STATUS=INEWI ) 

WRITE    (1,*)    NO_ TERMS 

DO   1=1  ,NO_  TERMS 

WRITE(*,300)    'expC-TAUCD/t)*' 
WRITE (lr*)   TAU(I),OMEGA(I) 
WRITE (1,*)   MULT(I) 

DO  J=1,MULT(I) 

WRITE(*,301)  '  C,RESP(I,0,J-1)  /f'^-l, 
+  'COS'  ,OMEGA(I),lt  +    ■ 

WRITE (1,*)   RESP(I,0,J-1) 

WRTTE(*,301)  '  '  ,RESP(I,1,J-1)  ,'t~',J-l, 
+  'SIN' ,OMEGA(I),'t)      ' 

WRITE (1,*)    RESP(I,1,J-1) 

ENDDO 
ENDDO 
CLOSE    (UNTr=l,    STATUS=IKEEP,) 

300  FORMAT      (A5,  E12.4E3/A3) 

301  FORMAT      (A2,E12.4E3  , A2, 12, A3, E12.4E3  ,A4) 

END 


******************************************************************* 

*  Department  of  Electrical  and  Computer  Engineering      * 

*  Kansas  State  University  * 

*  * 

*  VAX  FORTRAN  source  filename:     LOTS_0_ PLOTS. FDR  * 
***************************************************************** 

* 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


ROUTINE: 


DESCRIPTION: 


DOCUMENTATION 
FILES: 


ARGUMENTS: 
DEVICE 


subroutine 

FIRST,  PLOT  (DEVICE,  NUM_  POINTS,  X_DATA, 

Y_DATA,  X_AXIS_  TITLE,  X_AXIS_  UNITS, 

Y_AXIS_  TITLE,  Y_AXIS_  UNITS,  PLOT_TITLE, 

INFO) 


Makes  a  plot  using  X_DATA  as  abscissa 
and  Y_DATA  as  ordinate.     The  axes  are 
labelled  with  titles  and  units.     The 
plot  is  also  titled. 


None. 


(input)   integer 

is  the  device  type  to  display  the  plot 

7475  for  plotter 

4014  for  terminal   (Selanar  only) 


NUM_  POINTS     ( input )      intege  r 

is  the  number  of  data  points  to 
be  plotted 

X_DATA  (input)      real 

is  the  array  of  abscissa  values  for  the 
data  to  be  plotted 

Y_DATA  (input)     real 

is  the  array  of  ordinate  values  for  the 
data  to  be  plotted 

X_AXIS_TITLE   (input)     character* (*) 

is  the  title  to  be  placed  on  the  x_axis 

X_AXIS_UNITS   (input)     character* (*) 

is  the  name  to  be  given  to  the  units 


associated  with  the  x-axis 

Y_AXIS_TITLE   (input)     character* (*) 

is  the  title  to  be  placed  on  the  y-axis 

Y_AXIS_UNITS  (input)     character* (*) 

is  the  name  to  be  given  to  the  units 
associated  with  the  y-axis 

PLOT_  TITLE       (input)     character* (*) 

is  the  title  to  be  placed  on  the  plot 


RETURN: 
INFO 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

**************************************************************** 


(output)     real  (6) 

is  the  information  necessary  to  make 

subsequent  plots  on  the  same  axes. 


ROUTINES 
CALLED: 


P  System  of  Generalized  Plot  Routines 
AUTHOR:  James  F.    Stafford 

DATE  CREATED:        24May86     Version  1.0 


REVISIONS: 


None. 


SUBROUTINE     FIRST.  PLOT  (DEVICE ,  NUM_ POINTS ,  X_ DATA,  Y_DATA, 
+  X_  AXIS_  TITLE,  X_AXIS_  UNITS,  Y_AXIS_TITLEf 

+  Y_AXIS_ UNITS,  PLOT_TITLE,  INFO) 

IMPLICIT  NONE 

INTEGER     DEVICE,  NUM_POINTS,FORLAB,FORTIC,  NEGFLG,  FORM, 
+  SCNTL ,  LENSTR,  UPDOWN 

REAL  X_DATA(*)  ,Y_DATA(*)  , FACTOR, VEL,X,Y, LENGTH, 

+  FIRSTX,  DELTAX,  ANGLE,  CLEN,  FIRDEL  ( 4)  , 

+  DIVLNX,  DIVLNY,  WIDTH,  HEIGHT,  INFO  (6) 

CHARACTER*  (* )      X_ AXIS_TITLE,  Y_ AXIS_TITLE,  X_  AXIS_  UNITS, 
+  Y_  AXIS_  UNITS,  PLOT_ TITLE 


CHARACTER*  ( 1)      BLANK,  STZ  E 

*INTIALIZE  PLOT  DEVICE 

FACTOR=1.0 
BLANK='    ' 
SrZE='A' 

CALL  PINIT (DEVICE, BLANK, FACTOR,  SEE) 
*SET  PEN  VELOCITY 

VEL=10.0 

CALL     PSTVEL(VEL) 

♦ESTABLISH  ORIGIN 

X=4.5 
Y=4.5 

CALL  PORIG(X,Y) 

*SET  OFFSETS  FOR  AXIS  ROUTINES    (RELATIVE  TO  ORIGIN) 

X=O.0 
Y=0.0 

*DRAW  Y-AXIS  AND  LABEL 

LENGTH=12.0 

CALL  PSCALE(Y_DATA, NUM. POINTS,  LENGTH, FIRSTX, 
+  DELTAX,DIVLNY) 

FIRDEL(  3)  FIRSTX 
FIRDEL  ( 4)  =DELTAX 
FORLAB=110 
FORTIC=1001 
ANGLE=90.0 

CALL     PAXIS  (X,  Y,  Y_  AXIS_  TITLE,  Y_AXIS_  UNITS,  FORLAB, 
+  FORTIC,  LENGTH,  ANGLE,  FIRSTX,  DELTAX,  DIVLNY) 

*DRAW  X-AXIS  AND  LABEL 

LENGTH=18 

CALL  PSCALE(X_ DATA,  NUM. POINTS,  LENGTH, FIRSTX, 


+  DELTAX,DIVLNX) 

FIRDEL(1)=FIRSTX 

FIRDEL(2)=OELTAX 

F0RLAB=211 

F0RTIO2001 

ANGLE=0.0 

CALL  PAXIS  (X,  Y,  X_  AXIS_TITLE,  X_  AXIS_ UNITS,  POKLAB , 
+  FORTIC,  LENGTH,  ANGLE,  FIRSTX,  DELTAXf  DIVLNX) 

*DRiW  CURVE 

SCNTL=0 

CALL  PLINE(X_DATA,  Y_DATAf  NUM_  POINTS,  FIRDEL,  SCNTL, 
+  BLANK,  DIVLNX,DIVLNY) 

*TITLE  THE  PLOT 

UPDOWN=0 

X=9.0 

Y=13.0 

INFO(l)=FIRDEL(l) 

rNFO(2)=FIRDEL(2) 

INFO(3)=FIRDEL(3) 

INFO(4)=FIRDEL(4) 

INFO(5)=OIVLNX 

rNFO(6)=Dr7LNY 

CALL  PPLOT(X,Y,  UPTOWN) 

CALL  PTXTLN(PLOT_TITLE,LENSTR) 

WrDrH=-LENSTR/2 
HEIGHT=0.0 

CALL  PCHRPL (WIDTH,  HEIGHT) 

CALL  PTEXT (PLOT. TITLE) 

*  CALL  PCLOSP 

RETURN 

END 


***************************************************************** 

*  Department  of  Electrical  and  Computer  Engineering      * 

*  Kansas  State  University  * 

*  * 

*  VAX  FORTRAN  source  filename:         MULLER.FOR  * 
***************************************************************** 

* 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


ROUTINE: 


SUBROUTINE 
MULLER(POLYfDEG,X) 


DESCRIPTION:         This  program  uses  Muller'  s  method 
(page  262,  Numerical  Recipes)   to 
find  a  root  of  a  polynomial. 


DOCUMENTATION 

FILES :  None. 


ARGUMENTS: 


POLY 


DEG 


FACTOR 


NUM 


DEGF 


MULT 


RETURN: 


ROUTINES 
CALLED: 


The  following  arguments  are  passed 
to  the  subroutine: 

(input)   real 

is  an  array  containing  the  coefficients 

of  the  polynomial  of  interest. 

(input)    integer 

is  the  degree  of  the  above  polynomial. 

(input /out  put) 

is  an  array  containing  known  roots  of 

the  polynomial  represented  by  POLY. 

(input /out  put) 

is  the  number  of  factors  in  FACTOR. 

(input /out  put) 

is  an  array  containing  the  degree  of 

each  corresponding  factor  in  FACTOR. 

(input /out  put) 

is  an  array  containing  the  multiplicity 

of  each  corresponding  factor  in  FACTOR 

Not  used. 


DEFVAL,    COMPOSE 


* 
* 

*  AUTHOR:        James  F.  Stafford 

* 
* 

*  DATE  CREATED:        28May87  Version  1.0 
* 

* 

*  REVISIONS:  30Jun88  Added  deflation  and  factor  table 

*  updating. 
* 

* 
***************************************************************** 

SUBROUTINE  MULLER  (POLY,  DEC,  FACTOR,  NUM,  DEGF,  MULT) 

IMPLICIT  NONE 

INTEGER  DEG,  NUM,  DEGF  (*)  ,MULT(*)  ,  I,  NO_  ITERATIONS,  MAX 

REAL*8  POLY(0:*)  ,ZERO,FACTOR(10,0:2) 

COMPLEX*16  X(-2:l)  ,Q,A,B,C,D,P(-2:0)  ,DEFVAL 

PARAMETER  (ZERO=l  .OE-12) 

PARAMETER  (MAX=200) 

NO_  ITERATIONS=0 

X(-2)=DCMPLX(1.,1.) 
X(-1)=DCMPLX(1.,0.) 
X(0)=DCMPLX(1.,-1.) 

DO  WHILE    ((CDABS(X(0)-X(-1)).GT.CDABS(X(0))*2ERO) 
+  .AND.  (CDABS(X(0)-X(-2))  .GT.CDABS(X  (0)  )*ZERO) 

+  .  AND.  ( NQ_  ITERATIONS.  LT.  MAX) ) 

NO_  ITERATIONS=NO_  ITERATIONSH 

B=DCMPLX(0.,0.) 

D=DCMPLX(0.,0.) 

DO  WHILE    ((D.EQ.DCMPLX(0.,0.)).AND.  (B.  EQ.DCMPLX(0.  ,0.) ) ) 

DO  I=-2,0 

P (I )  =DEFVAL  (POLY,  DEG,  FACTOR,  NUM,  DEGF, MULT,  X  (I ) ) 

ENDDO 

Q=(X(0)-X(-l))/(X(-l)-X(-2)) 


A=Q*P(0)-Q*(l+Q)*P(-l)+Q*Q*P(-2) 

B=(2*C+l)*P(0)-((l+Q)**2)*P(-l)+Q*Q*P(-2) 

C=(1+Q)*P(0) 

D=SQRT(B*B-4*A*C) 

IF    ( (D. EQ. DCMPLX(0.  ,0.)). AND.  (B.EQ. DCMPLX(0.  ,0.)  ) )   THEN 

X(-l)=(X(0)+X(-l))/2. 
X(-2)=(X(0)+X(-2))/2. 

ENDIF 
ENDDO 
IF    (CDABS(Bt-D).GT.CDABS(B-D))   THEN 

X(l)=X(0)-(X(0)-X(-l))*2*C/(BfD) 
ELSE 

X(1)=X(0)-(X(0)-X(-1))*2*C/(B-D) 
ENDIF 
DO   I=-2,0 

X(I)=X(I+1) 

ENDDO 

ENDDO 

PRINT  VI  MADE   IT  HERE' /NO_  ITERATIONS 

NUM=NUMfl 

CALL  COMPOSE  (X  ( 1)  ,  FACTOR,  NUM,  DEGFf  MULT) 

RETURN 

END 


***************************************************************** 


DESCRIPTION: 


*  Department  °f  Electrical  and  Computer  Engineering  * 

*  Kansas  State  University  * 

*  * 

*  VAX  FORTRAN  source  filename:         PART.  FRAG  FOR  * 
************************************************************ 

* 

*  ROUTINE:  PROGRAM 

*  TEST 
* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

** 


DOCUMENTATION 
FILES: 


ARGUMENTS: 


This  program  computes  the  partial 
fraction  expansion  of  a  rational 
function  using  the  method  described  in 
the  thesis.     The  user  is  prompted 
for  a  filename  under  which  the  factored 
form  of  the  rational  function  has  been 
stored.     The  user  is  prompted  again 
for  a  filename  under  which  to  store  the 
partial  fraction  expansion. 


None. 
None. 


RETURN: 


ROUTINES 
CALLED: 


AUTHOR: 


Not  used. 


SPEC. READ,  PART. WRITE,  EXPAND 


James  F.    Stafford 


DATE  CREATED:       10Jun87  Version  1.0 


REVISIONS: 


None. 


************************************************************** 
PROGRAM  TEST 


IMPLICIT 


NONE 


INTEGER  I,  J,  K,DEGN,DEGD(10)  ,N0_ FACTS,  MULTS(IO)  , 

+  DEGX(10,5) 

REAL*8  NUM(0:15)  ,DEN(0:2,10)  ,X(0:1,10,5)  , 

+  FACT  (0:2) 

LOGICAL  EASY,  HARD 

CALL  SPEC_ READ  ( NUM,  DEGN,  DEN,  DEGD,  MULTS,  NO.  FACTS ) 

CALL  EXPAND ( NUM, DEGN,  DEN,  DEGD,  MOLTS,  X,DEGX,NO_ FACTS) 

DO  1=1,  NO_  FACTS 

DO  J=1,MULTS(I) 

PRINT  *,I,J 

DO  K=0,DEGX(I,J) 

PRINT  *,X(K,I,J) 

ENDDO 

ENDDO 

ENDDO 

CALL  PART_ WRITE  (NO_ FACTS,  MULTS,  X,  DEGX,  DEN,  DEGD) 

END 


***************************************************************** 

*  Department  of  Electrical  and  Computer  Engineering      * 

*  Kansas  State  University  * 

*  * 

*  VAX  FORTRAN  source  filename:         PART. WRITE.  FOR  * 
***************************************************************** 

* 

ROUTINE:  SUBROUTINE 

PART.  WRITE  (NO_  FACTS,  MULTS,  X,  DEGX,  DEN, 
DEGD) 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


DESCRIPTION: 


DOCUMENTATION 
FILES: 


ARGUMENTS: 


NO  FACTS 


MULTS 


X 


DEGX 


DEN 


This  program  writes  the  partial  fraction 
expansion  into  a  file.     The  user  is 
prompted  for  a  filename. 


None. 


The  following  arguments  are  passed  to 
the  subroutine: 

(input)    integer 

is  the  number  of  factors  in  the 

denominator  polynomial. 

(input)    integer 
is  an  array  containing  the 
multiplicities  of  each  factor  in 
the  denominator  polynomial. 

(input)   real 

is  a  three-dimensional  array.     X(If JfK) 
represents  the  Ith  coefficient  of  the 
numerator  of  the  (J,  K)th  term  in  the 
partial  fraction  expansion.     Namely, 
that  term  with  the  Jth  factor  of  DEN 
to  the  Kth  power  as  denominator. 

(input)  integer 

is  an  array.  DEGX  (I,  J)  represents  the 
degree  of  the  numerator  of  the  (I,J)th 
term  in  the  partial  fraction  expansion. 
See  the  description  of  X. 

(input)    real 

is  a  two-dimensional  array.     DEN(I,J) 
represents  the  coefficient  of  the  Ith 
power  of  x  in  the  Jth  factor  of  the 
denominator  polynomial. 


DEGD 


RETUEN: 


ROUTINES 
CALLED: 


AUTHOR: 

DATE  CREATED: 

REVISIONS: 


(input)    integer 

is  an  array.     DEGD  (I)   represents  the 
degree  of  the  Ith  factor  of  the 
denominator  polynomial. 

Not  used. 


None. 

James  F.    Stafford 

7Jun87     Version  1.0 

None. 


* 

* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
************************************************************ 

SUBROUTINE  PART.  WRITE  (NO_  FACTS,  MULTSf  X,  DEGXf  DEN,  DEGD) 


IMPLICIT 
INTEGER 


NONE 

I,  J,  K, DEGD (10)  ,NO_FACTS,MULTS(10)  , 
DEGX(10,5) 


REAL*8 
CHARACTER*  15 


200 


DEN(0:2,10)  ,X(0:1,10,5)  ,ALPHA,BETA,  A, B 

FILENAME 

PRINT  *, '  Enter  filename.  • 

READ    (*f200)    FILENAME 

FORMAT    (A15) 

OPEN   (UNIT=1,   FILE=FILENAME,    STATUS^  NEW' ) 

WRITE    (1,*)    NO_ FACTS 

DO  1=1  ,NO_  FACTS 

WRITE    (1,*)   DEGD(I) 

PRINT  *,' ORDER  =  ',DEGD(I) 

WRITE    (1,*)    MULTS(I) 

PRINT  *,  '  MULTIPLICITY  =I/MULTS(I) 


IF    (DEGD(I).EQ.l)    THEN 

WRITE    (1,*)   DEN(0,I) 
PRINT  *,' ALPHA  =',DEN(0,I) 

ELSE 

ALPHA=DEN(l,I)/2 
BETA=DSQRT(DEN(  0, 1)  -ALPHA*  *2) 
WRITE    (1,*)   ALPHA 
PRINT  V ALPHA  =', ALPHA 
WRITE    (1,*)    BETA 
PRINT  *, '  BETA  = ' ,BETA 

ENDIF 

DO  J=l,MOLTS(I) 

IF    (DEGD(I).EQ.l)    THEN 

WRITE    (1,*)    X(0,I,J) 
PRINT  *,  'A  =',X(0,I,J) 


ELSE 


A=X(1,I,J) 

B=  ( X  ( 0 , 1 ,  J)  -A*ALPH A)  /BETA 

WRITE    (1,*)   A 

PRINT  *,*A  ='fA 

WRITE    (1,*)    B 

PRINT  VB  =',B 


ENDIF 
ENDDO 
ENDDO 

CLOSE    (UNIT=1,    STATUS=IKEEP' 
RETURN 
END 


***************************************************************** 

*  Department  of  Electrical  and  Computer  Engineering      * 

*  Kansas  State  University  * 

*  * 

*  VAX  FORTRAN  source  filename:         PLOT. FOR  * 
***************************************************************** 

* 

*  ROUTINE:  PROGRAM 

PLOT 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

** 


DESCRIPTION: 


DOCUMENTATION 
FILES: 


ARGUMENTS: 


RETURN: 


ROUTINES 
CALLED: 


AUTHOR: 


This  program  makes  plots  of  time 
domain  response  functions  computed 
by  the  inverse  transform  program. 


None. 


None 


Not  used. 


FIRST. PLOT,  READ,  PLOT, 

PCLOSP  (contained  in  the  P  System 

of  Generalized  Plotting  Routines) 

James  F.    Stafford 


DATE  CREATED:        24May8  7  Version  1.0 


REVISIONS: 


None. 


************************************************************** 
PROGRAM  TEST 


IMHilCIT 

INTEGER 

REAL 


NONE 

DEVICE,  NUM_ POINTS,  I,  NUM_ FILES 

X_DATA(1000)  ,Y_DATA(1000)  , 
INFO  (6)  ,TONE,TTWO 


CHARACTER*  ( 15)      X_ TITLE,  Y_ TITLE,  X_  UNITS,  Y_  UNITS, 
+  TITLE,  FILES  (5) 


PRINT  *,'  INPUT  <7475>  FOR  PLOTTER  OR  <4014>  FOR  TERMINAL' 

READ    (*,*)DEVICE 

NUM_POINTS=1000 

X_TrrLE='TIME' 

Y_TITLE= 'VALUE' 

X_UNITS= '  INTERVALS' 

Y_UNTTS=' UNITS' 

TnLE='TEST  PLOT' 

PRINT  *,'  Input  number  of  files  to  plot' 
READ    (*,*)    NUM. FILES 

DO  1=1,  NUM_  FILES 

PRINT  *,'  Input  name  of  file  number  '  ,1 
READ    (*,200)    FILES (I) 


ENDDO 
200     FORMAT  (A15) 


PRINT  *,'  Input  initial  time' 
READ    (*,*)   TONE 

PRINT  *,'  Input  final  time' 
READ   (*,*)   TTWO 

DO  1=1  ,NUM_  POINTS 

X_DATA(I )  =TONE+  (TTWO-TONE)  * 

( FLOATJ  ( 1-1 )  /FLOATJ  ( NUM_  PO INTS ) ) 

ENDDO 

CALL  READ (NUM_ POINTS,  X_DATA,  Y_DATA,  FILES  ( 1) ) 

CALL  FIRST.  PLOT  (DEVICE,  NUM.  POINTS,  X_DATA,  Y.DATA, 
X_ TITLE,  X_ UNITS,  Y_TITLE,Y_  UNITS,  TITLE,  INFO) 

DO  1=2,  NUM_  FILES 

CALL  READ  ( NUM.  POINTS,  X_DATA,  Y_  DATA,  FILES  ( I ) ) 
CALL  PLOT  (DEV  ICE ,  NUM_  PO  INTS ,  X_  DATA,  Y_  DATA,  INFO ) 

ENDDO 

CALL  PCLQSP 


END 


***************************************************************** 

*  Department  of  Electrical  and  Computer  Engineering  * 

*  Kansas  State  University  * 

*  * 

*  VAX  FORTRAN  source  filename:     PLOT. Q_MATIC. FOR  * 
***************************************************************** 

* 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


ROUTINE: 


DESCRIPTION: 


DOCUMENTATION 
FILES: 


ARGUMENTS: 
DEVICE 


subroutine 

PLOT  (DEV  ICE ,  NUM_  PO INTS ,  X_  DATA, 

Y_DATA,  INFO) 


Makes  a  plot  using  X_DATA  as  abscissa 
and  Y_DATA  as  ordinate.     The  axes  are 
assumed  to  be  already  drawn  in 
accordance  with  INFO. 


None. 


(input)   integer 

is  the  device  type  to  display  the  plot 

7475  for  plotter 

4014  for  terminal   (Selanar  only) 


NUM_POINTS     (input)      integer 

is  the  number  of  data  points  to 
be  plotted 


X  DATA 


Y  DATA 


RETURN: 
INFO 


(input)     real 

is  the  array  of  abscissa  values  for  the 

data  to  be  plotted 

(input)     real 

is  the  array  of  ordinate  values  for  the 

data  to  be  plotted 


(output)     real  (6) 

is  the  information  necessary  to  make 

subsequent  plots  on  the  same  axes. 


* 

*  ROUTINES 

*  CALLED:  P  System  of  Generalized  Plot  Routines 

* 

* 

*  AUTHOR:        James  F.  Stafford 

* 

*  DATE  CREATED:        24May86     Version  1.0 
* 

* 

*  REVISIONS:  None. 
* 

* 
**************************************************************** 


SUBROUTINE     PLGT(DEVICE,  NUM_  POINTS,  X_DATA,  Y_DATA, 
+  INFO) 

IMPLICIT  NONE 

INTEGER     DEVICE,  NUM_  POINTS,  FORLABrFORTICrNEGFLGf  FORM, 
+  SCNTL ,  LENSTR,  UPDOWN 

REAL  X_DATA(*)  ,Y_DATA(*)  ,FACTOR,VEL,X,Y, LENGTH, 

+  FIRSTX,DELTAX,ANGLE,CLEN,FIRDEL(4)  , 

+  DIVLNX,DIVLNY,  WIDTH,  HEIGHT,  INFO  (6) 

CHARACTER*  ( 1)      BLANK,  SE  E 

*INTIALIZE  PLOT  DEVICE 

FACTOR=1.0 
BLANK='    ' 
SIZF^'A' 

*  CALL     PINIT (DEVICE, BLANK, FACTOR,  SEE) 

*SET  PEN  VELOCITY 

VEL=10.0 

CALL     PSWEL(VEL) 
♦ESTABLISH  ORIGIN 

X=4.5 

Y=4.5 


CALL  PORIG(X,Y) 

*SET  OFFSETS  FOR  AXIS   ROUTINES    (RELATIVE  TO  ORIGIN) 

X=0.0 
Y=0.0 

♦ESTABLISH  INFORMATION  FOR  DOTTING  SUBROUTINE 

FIRDEL(1)=INF0(1) 

FIRDEL(2)=INFO(2) 

FIRDEL(3)=INFO(3) 

FIRDEL(4)=INFO(4) 

DIVLNX=INFO(5) 

DIVLNY=INFO(6) 

*DRAW  CURVE 

SCNTL=0 

CALL  PL INE  (X_DATA,  Y_  DATA,  NUM_  PO INTS ,  FIRDEL ,  SCNTL , 
+  BLANK,  DIVLNX,DIVLNY) 

RETURN 

END 


***************************************************************** 


* 

* 

* 

* 

** 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


Department  of  Electrical  and  Computer  Engineering 
Kansas  State  University 


* 
* 
* 

VAX  FORTRAN  source  filename:         POLADD.FOR  * 

*************************************************************** 


ROUTINE: 


DESCRIPTION: 


DOCUMENTATION 
FILES: 


ARGUMENTS: 


DEGA 


B 


DEGB 


RETOBN: 


ROUTINES 
CALLED: 


SUBROUTINE 

POLADD  (A,  DEGA,  B,  DEGB) 


This  program  subtracts  on  polynomial 
from  another,    replacing  the  first  addend 
with  the  difference. 


None. 


The  following  arguments  are  passed  to 
the  subroutine. 

( input /out  put )   r  eal 

is  an  array  containing  the  coefficients 

for  a  polynomial  on  input  and  the 

coefficients  for  the  difference  on 

return. 

( input/output )    integer 
is  the  degree  of  the  above  polynomial 
on  input  and  the  degree  of  the 
difference  on  return. 

(input)   real 

is  an  array  containing  the  coefficients 
of  the  polynomial  to  be  subtracted 
from  A. 

(input)  integer 

is  the  degree  of  the  above  polynomial. 


Not  used. 


None. 


*  AUTHOR:        James  F.  Stafford 
* 

* 

*  DATE  CREATED:        6Jun87     Version  1.0 

* 
* 

*  REVISIONS:  None. 
* 

* 
**************************************************************** 

SUBROUTINE  FOLADD  (A,  DEGAf  B,  DEGB) 

IMPLICIT  NONE 

INTEGER  I ,  DEGA,  DEGB ,  ADDS 

REAL*8  A(0:*),B(0:*) 

ADDS=0 

IF    (DEGA.  GE.  DEGB)    THEN 

DO   1=0 ,  DEGB 

A(I)=A(I)-B(I) 
ADDS=ADDS+1 

ENDDO 

ELSE 

DO  1=0,  DEGA 

A(I)=A(I)-B(I) 
ADDS=ADDS+1 

ENDDO 

DO  I=DEGA+1,DEGB 

A(I)=-B(I) 
ADDS=ADDS+1 

ENDDO 

ENDIF 

DEGA=JMAX0  (DEGA,  DEGB) 

PRINT  *, ADDS, 'additions' 


RETURN 
END 


***************************************************************** 


Department  of  Electrical  and  Computer  Engineering 
Kansas  State  University 


DESCRIPTION: 


DOCUMENTATION 
FILES: 


ARGUMENTS: 


* 
* 

*  * 

*  VAX  FORTRAN  source  filename:         POLDIV.FOR  * 
***************************************************************** 

* 

*  ROUTINE:  SUBROUTINE 

*  POLDIV(NUMfN_DEGfDENfD_DEGfQUOfDEGQ, 

*  REM,  DEGRf  EASY,  HARD) 
* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


NUM 


N  DEG 


DEN 


D  DEG 


QUO 


DEGQ 


REM 


This  program  performs  the  division 
algorithm  on  two  input  polynomials 
to  obtain  a  quotient  and  remainder. 


None. 


The  following  arguments  are  passed  to 
the  subroutine. 

(input)  real 

is  an  array  containing  the  coefficients 

of  the  numerator  polynomial. 

(input)    integer 

is  the  degree  of  the  numerator  polynomial. 

(input)   real 

is  an  array  containing  the  coefficients 

of  the  denominator  polynomial. 

(input)    integer 

is  the  degree  of  the  denominator 

polynomial. 

(output)   real 

is  an  array  containing  the  coefficients 

of  the  quotient  polynomial. 


nomial, 


(output)    integer 

is  the  degree  of  the  quotient  poly 

(output)    real 

is  an  array  containing  the  coefficients 

of  the  remainder  polynomial. 


* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 

* 
* 
* 
* 
* 
* 


DEGR 


EASY 


HARD 


RETURN: 


ROUTINES 
CALLED: 


AUTHOR: 


(output)    integer 

is  the  degree  of  the  remainder 

polynomial. 

(input)  logical 

should  be  set  to  .TRUE,  if  both  the 
numerator  and  denominator  are  monic 
polynomials.  This  will  save  calculations. 

(input)  logical 

must  be  set  to  .TRUE,  if  the 

denominator  is  not  a  monic  polynomial. 

Otherwise  leave  it  false  to  save 

calculations. 


Not  used. 


None. 


James  F.   Stafford 


DATE  CREATED:        6Jun87     Version  1.0 


REVISIONS: 


27Jul87  Added  calculation-saving. 


**************************************************************** 

SUBROUTINE  POLDIV  ( NUM,  N_DEG,  DEN,  D_DEGf  QUO,  DEGQ, 

+  REM,  DEGR,  EASY,  HARD) 


IMPLICIT 


NONE 


INTEGER  J,K,N_DEG,D_DEE,DB^,DEGR,MULT,ADD,DIV 

REAL*8  NUM(0:*) ,DEN(0:*) ,QUO(0:*) ,REM(0:*) ,ZERO 

LOGICAL  EASY, HARD 

PARAMETER  ( Z  ERO=l . 0  E-5 ) 

DEGQ=N_DEG-D_DEG 

MULT=0 

DIV=0 

ADD=0 


DO  J=0,N_DEG 

REM(J)=NUM(J) 
ENDDO 
IF    (DEGQ.GE.O)    THEN 

DO  K=DEGQ,0,-1 

IF    (HARD)   THEN 

QUO(K)  =REM(D_DEGf  K)/DEN(D_DEG) 
DIV=DIV+1 

ELSE 

QUO(K)=REM(D_DEG+K) 
ENDIF 
DO  J=D_DEG+K-1,K,-1 

IF    (EASY)   THEN 

REM(J)  =REM(J)  -DEN(  J-K) 
ADD=ADDfl 

ELSE 

REM  (J)  =REM  (J)  -QUO  (K )  *DEN  ( J-K ) 

ADD=ADDH 

MULT=MULT+1 

ENDIF 

ENDDO 

EASY=.  FALSE. 
ENDDO 

IF    (D_DEG.EQ.O)    REM(0)=0. 
DEGR=  JMAXO  ( 0 ,  D_  DEG-1 ) 
DO  WHILE    ( (DAB  S  ( REM  (DEGR)  )  .LT.  ZERO).  AND.  DEGR.CT.0) 

DEGR=DEGR-1 


ENDDO 

ELSE 

DEGR=N_DEG 

DEGQ=0 

QUO(0)=0. 

ENDIF 

*  PRINT  *,njLT, 'multiplies1 

*  PRINT  * , ADD, '  addi  t i  ons ' 

*  PRINT  *,DIV, 'divisions' 

RETURN 
END 


***************************************************************** 

*  Department  of  Electrical  and  Computer  Engineering  * 

*  Kansas  State  University  * 

*  * 

*  VAX  FORTRAN  source  filename:         FOLMULT.  FOR  * 
***************************************************************** 

* 

ROUTINE:  SUBROUTINE 

FOLMULT  (FOL1 ,  DEGl ,  FOL2 ,  DEG2 ,  PROD,  DEGP, 
EASY) 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


DESCRIPTION: 


DOCUMENTATION 
FILES: 


ARGUMENTS: 


P0L1,P0L2 


DEG1,DEG2 


PROD 


DEGP 


EASY 


RETURN: 


ROUTINES 
CALLED: 


This  program  multiplies  two  polynomials 
and  returns  their  product. 


None. 


The  following  arguments  are  passed  to 
the  subroutine. 

(input)   real 

are  arrays  containing  the  coefficients 

of  the  polynomials  to  be  multiplied. 

(input)    integer 

are  the  degrees  of  the  above 

polynomials. 

(output)   real 

is  an  array  containing  the  coefficients 

of  the  product  polynomial. 

(output)    integer 

is  the  degree  of  the  product 

polynomial. 

(input)    logical 

should  be  set  to  .TRUE,   only  if  POL1 

is  monic  in  order  to  save  calculations. 


Not  used. 


None. 


*  AU1H0R:        James  F.  Stafford 
* 

* 

*  DATE  CREATED:   7Jun37  Version  1.0 
* 

* 

*  REVISIONS:  20Jul87  Added  calculation- saving. 
* 

* 

**************************************************************** 

SUBROUTINE  POLMULT(POLl , DEGl ,  FOL2 , DEG2 ,  PROD,  DEGP,  EASY) 

IMPLICIT  NONE 

INTEGER  J,  K, DEGl, DEG2, DEGP,  ADD,  MULT 

REAL*8  FOLl(0:*)fPOL2(0:*),PRCD(0:*) 

LOGICAL  EASY 

ADD=0 

MULT=0 

DEGP=DEG1+DEG2 

DO  J=DEG2,0,-1 

IF    (EASY)   THEN 

PRCD(DEG1+J)  =POL2  ( J) 

ELSE 

PRCD(DEG1+J)  =P0L1  (DEGl)  *POL2  ( J) 
MULT=MULT+1 


ENDIF 
ENDDO 
DO  J=DEG1-1,0,-1 

DO  K=€)EG2,1,-1 


PROD ( J+ K)  =PROD  ( J+ K) +POL1  ( J)  *POL2  (K) 

ADD=ADDfl 

MULT=MULT+1 


ENDDO 


FRCD  (J)  =P0L1  ( J)  *P0L2  ( 0) 
MULT=MULT+1 

ENDDO 

PRINT  *,  ADD,  ■  additions1 
PRINT  *, MULT,  'multiplies' 

RETURN 

END 


***************************************************************** 

*  Department  of  Electrical  and  Computer  Engineering      * 

*  Kansas  State  University  * 

*  * 

*  VAX  FORTRAN  source  filename:         POL Y_ READ. FOR  * 
***************************************************************** 

* 

* 

* 
* 
* 

* 
* 
* 

* 
* 
* 
* 

* 


ROUTINE: 


DESCRIPTION: 


DOCUMENTATION 
FILES: 


SUBROUTINE 

POLY.  READ  (POLY,  DEC ) 


This  program  reads  a  polynomial  from 
a  data  file.     The  user  must  enter 
the  data  file  name  from  the  keyboard. 


None. 


* 
* 
* 
* 

* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 


ARGUMENTS: 
POLY 

DEG 

RETURN: 


ROUTINES 
CALLED: 


AUTHOR: 


The  following  arguments  are  passed  to 
the  subroutine. 

(output)    real 

is  an  array  to  receive  the  coefficients 

of  the  polynomial. 

(output)    integer 

is  the  degree  of  the  polynomial. 


Not  used. 


None. 


James  F.   Stafford 


DATE  CREATED:        6Jun87     Version  1.0 


* 
* 
* 
* 
* 
* 
* 
**************************************************************** 


REVISIONS: 


None. 


SUBROUTINE 


POLY_  READ  (POLY,  DEG ) 


IMPLICIT 

NCNE 

INTEGER 

DBG,  I 

REAL*8 

POLY(0:10) 

READ    (1,*) 

DEG 

READ    (1,*) 

(POLY  (I),    1=0,  DEG) 

RETURN 

END 

***************************************************************** 

*  Department  of  Electrical  and  Computer  Engineering      * 

*  Kansas  State  University  * 

*  * 

*  VAX  FORTRAN  source  filename:         READ. FOR  * 
***************************************************************** 

* 

SUBROUTINE 

READ  (NUM. POINTS,  X_DATA,  Y_DATA,  FILENAME) 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


ROUTINE: 


DESCRIPTION: 


DOCUMENTATION 
FILES: 


ARGUMENTS: 


This  program  evaluates  a  response 
function  computed  by  the  inverse 
transform  program.     The  particular 
function  parameters  are  stored  in  a 
data  file  under  FILENAME. 


None. 


The  following  arguments  are  passed  to 
the  subroutine: 


NUM. PO INTS     ( input )    int ege r 

is  the  number  of  points  in  X_DATA 
at  which  response  values  are  desired. 


X  DATA 


Y  DATA 


FILENAME 


RETURN: 


ROUTINES 
CALLED: 


AUTHOR: 


(input)    real 

is  an  array  containing  the 
values  of  time  that  the  response 
function  is  to  be  evaluated  at. 

(output)   real 

is  an  array  to  accumulate  the 

computed  function  values. 

( input )   character 

is  the  filename  of  the  response 

function  to  be  evaluated. 


Not  used. 


None. 


James  F.    Stafford 


* 

* 

*  DATE  CREATED:        24May86  Version  1.0 

* 

*  REVISIONS:  None. 
* 

* 
**************************************************************** 

SUBROUTINE  READ  ( NUM_  POINTS,  X.DATA,  Y_DATAf  FILENAME ) 

IMPLICIT  NONE 

INTEGER  NUM_ POINTS,  I,  J,K,NO_ TERMS,  MULTS 

REAL  X_DATA(1000)  ,Y_DATA(1000)  , EXPONENTIAL, 

+  TONE,  TTWO,  POWEROFT 

REAL*8  TAU,  OMEGA,  RESP  (0:1) 

CHARACTER*  (15)      FILENAME 

DO  1=1,  NUM_  POINTS 

Y_DATA(I)=0. 

ENDDO 

OPEN   ( UNTT=1  ,FILE=FILENAME,  STATUS=  'OLD' ) 

READ    (1,*)    NO_ TERMS 

DO  1=1,  NQ_  TERMS 

READ(1,*)   TAU, OMEGA 
READ(1,*)   MULTS 

POWEROFT=l 

DO  J=l, MULTS 

READ(1,*)    RESP(O) 
READ(1,*)    RESP(l) 

DO  K=1,NUM_  POINTS 

IF    (J-1.GT.0)    POWEROFT=X_DATA(K)**(J-1) 

Y_DATA  (K )  =Y_  DATA  (K )  + 
+  SN3L(DEXP(-TAU*X_DATA(K))* 


+  (RESP(0)*DCDS(OMEGA*X_DATA(K))  + 

+  RESP(1)*DSIN(0MBGA*X_DATA(K))))* 

+  POWEROFT 

ENDDO 

ENDDO 

EM)DO 

CLOSE   (UNIT=1) 

RETURN 

END 


***************************************************************** 

*  Department  of  Electrical  and  Computer  Engineering      * 

*  Kansas  State  University  * 

*  * 

*  VAX  FORTRAN  source  filename:         ROOT. FIND.  FDR  * 
***************************************************************** 

* 

PROGRAM 
TEST 


ROUTINE: 


DESCRIPTION: 


DOCUMENTATION 
FILES: 


ARGUMENTS: 


This  program  factors  the  denominator  of 
a  rational  function  with  real  coefficients 
into  irreducible  polynomials  in  R[x]. 
The  user  is  prompted  for  a  filename 
under  which  a  rational  function  has 
been  stored.     The  user  is  again  prompted 
for  a  filename  to  store  the  result 
under. 


None. 


* 
* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

**************************************************************** 

PROGRAM  TEST 


RETURN: 


ROUTINES 
CALLED: 


AUTHOR: 


None. 


Not  used. 


POL Y_ READ,  FACTORER,  SPEC. WRITE 


James  F.   Stafford 


DATE  CREATED:       8Jun37     Version  1.0 


REVISIONS: 


None. 


IMPLICIT 

INTEGER 

REAL*8 


NONE 

DEGN,DEGD,DEGF(10)  ,MJM_  FACTS,  MULT  (10)  ,1,  J 

NUM(0:10)  ,DENOM(0:10)  ,FACTOR(10,0:2) 


CHARACTER* 15         FILENAME 
PRINT  *,'  Input  data  file  name' 
READ    (*,200)    FILENAME 
200  FORMAT    (A15) 

OPEN    ( UNIT=1 ,  FILE=FILENAMEf  STATUS=  *  OLD' ) 

CALL  POL Y_  READ  (NUM,DEGN) 

CALL  POL Y_ READ  (DENOMfDEGD) 

CLOSE    (UNTT=1) 

CALL  FACTORER  (DENOMf  DEGD,  FACTOR,  NUM_  FACTS,  DEGF,  MULT) 

DO  1=1  ,NUM_  FACTS 

PRINT  V  FACTOR  NUMBER  '  rI 

DO  J=0,DEGF(I) 

PRINT  *,  FACTOR  (I,  J) 

ENDDO 

PRINT  ♦/MULTIPLICITY1  ,  MULT  (I) 
ENDDO 

CALL  SPEC. WRITE  ( NUM,  DEGN,  NUM_ FACTS,  FACTOR, DEGF,  MULT) 
END 


***************************************************************** 

*  Department  of  Electrical  and  Computer  Engineering      * 

*  Kansas  State  University  * 

*  * 

*  VAX  FORTRAN  source  filename:     SIMFLE._ILOT.FOR  * 
***************************************************************** 


* 
* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


ROUTINE: 


DESCRIPTION: 


DOCUMENTATION 
FILES: 


ARGUMENTS: 
DEVICE 


subroutine 

SIMPLE.  PLOT  (DEVICE,  NUM_  POINTS,  X_DATA, 
Y_  DATA,  X_  AXIS_  TITLE,  X_  AXIS_  UNITS, 
Y_AXIS_ TITLE,  Y_AXIS_  UNITS,  PLOT. TITLE, 
PLOT.  TYPE) 


Makes  a  plot  using  X_DATA  as  abscissa 
and  Y_DATA  as  ordinate.     The  axes  are 
labelled  with  titles  and  units.     The 
plot  is  also  titled.     One  of  four  plot 
types  can  be  selected. 


None. 


(input)   integer 

is  the  device  type  to  display  the  plot 

7475  for  plotter 

4014  for  terminal    (Selanar  only) 


NUM_  POINTS     (input)      integer 

is  the  number  of  data  points  to 
be  plotted 

X_DATA  (input)      real 

is  the  array  of  abscissa  values  for  the 
data  to  be  plotted 

Y_DATA  (input)      real 

is  the  array  of  ordinate  values  for  the 
data  to  be  plotted 

X_AXIS_TITLE  (input)     character* (*) 

is  the  title  to  be  placed  on  the  x_axis 

X_AXIS_ UNITS   (input)     character* (*) 


is  the  name  to  be  given  to  the  units 
associated  with  the  x-axis 

Y_AXIS_  TITLE  (input)     character* (*) 

is  the  title  to  be  placed  on  the  y-axis 

Y_AXIS_ UNITS  (input)     character* (*) 

is  the  name  to  be  given  to  the  units 
associated  with  the  y-axis 


PLOT  TITLE 


PLOT  TYPE 


( input )     character*  (* ) 

is  the  title  to  be  placed  on  the  plot 

(input)     character* (*) 
is  a  character  string  which 
specifies  the  type  of  plot  to  be 
generated.     The  following  are  valid: 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

**************************************************************** 


RETURN: 


ROUTINES 
CALLED: 


AUTHOR: 


'LINEAR' 
'LOG_  LINEAR1 
'LINEAR.  LOG' 
'LOGLLOG' 

Not  used. 


for  linear- linear 
for  log-linear 
for  linear-log 
for  log-log 


P  System  of  Generalized  Plot  Routines 


James  F.   Stafford 


DATE  CREATED:        5Sep36     Version  1.0 


RB/ IS  IONS: 


None. 


SUBROUTINE     SIMPLE. PLOT  (DEVICE,  NUM.  POINTS,  X_ DATA,  Y.DATA, 
+  X_  AXIS_  TITLE,  X_AXIS_  UNITS,  Y_AXIS_  TITLE, 

+  Y_  AXIS_  UNITS,  PLOT. TITLE,  PLOT_ TYPE ) 


IMPLICIT  NONE 

INTEGER     DEVICE,  NUM_POINTS,FORLAB,PORTIC,NEGFLG,  FORM, 


+  SCNTL,LENSTR,  UPDOWN 

REAL  X_DATA(*)  ,Y_DATA(*)  ,  FACTOR,  VEL,X,  Y,  LENGTH, 

+  FIRSTX,  DELTAX,  DIVLEN,  ANGLE,  CLEN,  FIRDEL  ( 4)  , 

+  DIVLNX,  DTVLNY,  WIDTH,  HEIGHT 

CHARACTER*  (*)      X_AXIS_ TITLE,  Y_  AXIS. TITLE,  X_AXIS_  UNITS, 
+  Y_  AXIS_  UNITS ,  PLOT. TITLE,  PLOT.  TYPE 

CHARACTER*  (1)      BLANK, SIZE 

*INTIALIZE  PLOT  DEVICE 

FACTOR=1.0 
BLANK='    ' 
SIZF^'A1 

CALL     PINIT  (DEVICE,  BLANK,  FACTOR,  SJZ  E) 
*SET  PEN  VELOCITY 

VEL=10.0 

CALL     PSTVEL(VEL) 

♦ESTABLISH  ORIGIN 

X=4.5 
Y=4.5 

CALL  PORIG(X,Y) 

*SET  OFFSETS  FOR  AXIS  ROUTINES    (RELATIVE  TO  ORIGIN) 

X=0.0 
Y=0.0 

*DRAW  Y-AXIS  AND  LABEL 

IF    (PLOT. TYPE.  EQ.  ■  LINEAR'  .OR.  PLOT. TYPE.  EQ. '  LINEAR, LOG1  ) 
+  THEN 

LENGTH=12.0 

CALL  PSCALE  (Y_DATA,  NUM.  POINTS,  LENGTH,  FIRSTX, 
+  DELTAX,  DIVLEN) 

FIRDEL  ( 3)  =FIRSTX 
FIRDEL  ( 4)  =DELTAX 
DTVLNY=DIVLEN 


FORLAB=110 

PORTIOlOOl 

ANGLE=90.0 

CALL     PAXIS(X,Y,Y_AXIS_  TITLE,  Y_AXIS_  UNITS,  FORLAB, 
+  FORTIC,  LENGTH,  ANGLE,  FIRSTX,  DELTAX,  DIVLEN) 

ELSE 

LENGTH=12.0 

CALL  PLOGSC(Y_DATA,  NUM. POINTS,  LENGTH,  FIRSTX,CLEN, 
+  NEGFLG) 

FIRDEL(3)=FIRSTX 
FIRDEL(4)=CLEN 
FORM=-1010 
ANGLE=90.0 

CALL  PLGAXS(X,Y,Y_AXIS_TITLE,Y_AXIS_  UNITS,  FORM, 
+  LENGTH,  ANGLE,  FIRSTX,  CLEN) 

ENDIF 

*DRAW  X-AXIS  AND  LABEL 

IF    (PLOT. TYPE.  EQ.  ■  LINEAR1  .OR.  PLOT. TYPE.  EQ. '  LOCLINEAR' ) 
+     THEN 

LENGTH=18 

CALL  PSCALE  (X_  DATA,  NUM_  PO INTS,  LENGTH ,  FIRSTX, 
+  DELTAX,DIVLEN) 

FIRDEL(1)=FIRSTX 

FIRDEL(2)=€»ELTAX 

DIVLNX=DIVLEN 

FORLAB=211 

FORTI02001 

ANGLE=0.0 

CALL  PAXIS(X,Y,X_AXIS_ TITLE,  X_AXIS_ UNITS,  PORLAB, 
+  FORTIC,  LENGTH ,  ANG  LE,  FIRSTX,  DELTAX,  DIVLEN) 

ELSE 

LENGTH=18.0 

CALL  PLOGSC  (X_ DATA,  NUM_ POINTS,  LENGTH,  FIRSTX,  CLEN, 
+  NEGFLG) 


FIRDEL  ( 1)  =FIRSTX 
FIRDEL  ( 2)  =CLEN 
FORM=+2011 
ANGLE=0.0 

CALL  PLGAXS(X,Y,X_mS_  TITLE,  X_AXIS_  UNITS,  FORM, 
+  LENGTH,  ANGLE,  FIRSTX,CLEN) 

ENDIF 

*DFJW  CUR/E 

IF    (PLOT_TYPE.EQ.'  LINEAR') THEN 

SCNTL=0 

CALL  PLINE  (X_DATA,  Y_DATA,  NQM_ POINTS,  FIRDEL,  SCNTL, 
+  BLANK,  DIVLNX,DIVLNY) 

ELSE  IF  (PLOT.  TYPE.  EQ.'LOG_  LINEAR'  )THEN 

SCNTL=0 

CALL  PLGLIN  (X_DATA,  Y.DATA,  NUM_ POINTS,  FIRDEL,  SCNTL, 
+  BLANK,  DIVLEN) 

ELSE  IF  (PLOT. TYPE.  EQ.1  LINEAR. LOG1) THEN 

SCNTL=0 

CALL  PLNLCG  (X_DATA,  Y_DATA,  NUM. POINTS,  FIRDEL,  SCNTL, 
+  BLANK,  DIVLEN) 

ELSE  IF  (PLOT.  TYPE.  EQ.'LCG_  LOG')  THEN 

SCNTL=0 

CALL  PLGLOG(X_ DATA,  Y_ DATA,  NUM_ POINTS,  FIRDEL,  SCNTL, 
+  BLANK) 

ENDIF 

*TITLE  THE  PLOT 

UPDOWN=0 

X=9.0 

Y=13.0 

CALL  PPLOT(X,Y,UPDOWN) 


CALL  PTXTLN( PLOT. TITLE,  LENSIR) 

WIDTH=-LENSTR/2 
HEIGHT=0.0 

CALL  PCHRPL  (WIDTH,  HEIGHT) 

CALL  PTEXT  (PLOT. TITLE) 

CALL  paosp 

RETURN 
END 


***************************************************************** 


*  Department  of  Electrical  and  Computer  Engineering  * 

*  Kansas  State  University  * 

*  * 

*  VAX  FORTRAN  source  filename:  SPEC. INPUT. FDR  * 
***************************************************************** 

* 

PROGRAM 

SPEC  INPUT 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

• 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


ROUTINE: 


DESCRIPTION: 


DOCUMENTATION 
FILES: 


ARGUMENTS: 


RETURN: 


ROUTINES 
CALLED: 


AUTHOR: 


DATE  CREATED: 


REVISIONS: 


This  program  allows  the  user  to  enter 
the  specifications  for  a  rational 
function  already  in  factored  form. 


None. 
None. 
Not  used. 

None. 

James  F.    Stafford 

15Jun87  Version  1.0 

None. 


**************************************************************** 

NONE 


IMPLICIT 
INTEGER 

REAL 


I,  J,  N_DEG,D_DEG(10)  ,NO_RCOTS,MULTS(10)  , 
DEG,NO_  FACTS 

NUM(0:15)  ,DEN(0:2,10)  ,X(0:1,10,5)  , 
FACT (0:2)  ,FOLY(0:20) 


PRINT  *,  ■  Enter  numerator  information. ' 
CALL   INPUT. NCNFACT(NUM,N_DEG) 


PRINT  *, '  Enter  denominator  inf  ormation. ' 

PRINT  *,'  Input  number  of  relatively  prime  irreducible  factors.1 

READ   (*,*)    NO. FACTS 

DO  1=1,  N0_  FACTS 

PRINT  *,'  Enter  information  on  factor  number'  ,1 

CALL  INPUT.  NCNFACT  (FACT,  D_  DBG  (I)) 

DO  J=0,D_DEG(I) 

DEN(JfI)=FACT(J) 

ENDDO 

PRINT  *, ' Enter  number  of  times  this  factor  appears. ' 

READ  (*,*)  MULTS(I) 

ENDDO 

CALL  SPEC. WRITE  ( NUM,  N_ DEG ,  N0_ FACTS,  DEN,  D_ DEG ,  MULTS ) 

END 

SUBROUTINE  INPUT.  NCNFACT  (POLY,  DEG) 

IMPLICIT  NONE 

INTEGER  DEG,  I 

REAL  POLY(0:*) 

PRINT  *, '  Input  degree1 

READ    (*,*)   DEG 

DO  1=0,  DEG 

PRINT  *, 'Input  coeff.   of  power  ',1 
READ   (*,*)    POLY(I) 

ENDDO 

RETURN 

END 


DESCRIPTION: 


DOCUMENTATION 
FILES: 


ARGUMENTS: 


SUBROUTINE 

SPEC_READ  (NUM,  DEGN,  DEN,  DEGDf  MULTS, 

NO_  FACTS) 


This  program  reads  data  for  a  rational 
function  out  of  a  file  created  by  the 
program  for  factoring  the  denominator. 
There  is  also  a  program  called 
SPEC. WRITE  that  will  create  a  data  file 
in  the  proper  format. 


None. 


***************************************************************** 

*  Department  of  Electrical  and  Computer  Engineering  * 

*  Kansas  State  University  * 

*  * 

*  VAX  FORTRAN  source  filename:         SPEC_ READ. FOR  * 
***************************************************************** 

* 

*  ROUTINE: 
* 
* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

*  DEGD 
* 

* 

* 

* 

*  MULTS 
* 

* 


NUM 


DEGN 


DEN 


The  following  arguments  are  passed  to 
the  subroutine: 

(output)   real 

is  an  array  containing  the  coefficients 

of  the  numerator  polynomial. 

(output)  integer 

is  the  degree  of  the  numerator 

polynomial. 

(output)   real 

is  a  two-dimensional  array.     DEN(If  J) 
represents  the  coefficient  of  the  Ith 
power  of  x  in  the  Jth  factor  of  the 
denominator  polynomial. 

(output)    integer 

is  an  array.     DEGD(I)   represents  the 
degree  of  the  Ith  factor  in  the 
denominator  polynomial. 

(output)    integer 

is  an  array.     MULTS  (I)    represents  the 

multiplicity  of  the  Ith  factor  in  the 


NO  FACTS 


RETURN: 


ROUTINES 
CALLED: 


AUTHOR: 


DATE  CREATED: 


REVISIONS: 


denominator  polynomial. 

(output)  integer 

is  the  number  of  factors  in  the 

denominator  polynomial. 


Not  used. 

None. 

James  F.   Stafford 

6Jun87     Version  1.0 

None. 


* 
* 
* 

* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
**************************************************************** 

SUBROUTINE  SPEC. READ  ( NUM,  DEGN, DENf  DEGD,  MULTS,  NO_ FACTS ) 


IMPLICIT 

INTEGER 

REAL*8 


NONE 

I,  J,  DEGN,  DEGD (10)  ,NO_ FACTS,  MULTS  (10) 

NUM(0:15)  ,DEN(0:2,10)  ,X(0:1,10,5)  , 
FACT  (0:2) 


200 


CHARACTER*15         FILENAME 

PRINT  *, '  Input  data  file  name' 

READ   (*,200)    FILENAME 

FORMAT    (A15) 

OPEN   ( UNTT=1 ,  FILE=FILENAME,  STATUS= '  OLD1 ) 

READ    (1,*)   DEGN 

DO   1=0, DEGN 

READ    (1,*)    NUM(I) 
ENDDO 


READ  (1,*)  NO. FACTS 
DO  1=1,  NO_  FACTS 

READ  (1,*)  DEGD(I) 

DO  J=0,DEGD(I) 

READ  (1,*)  DEN (J, I) 

ENDDO 

READ  (1,*)  MULTS(I) 
ENDDO 
CLOSE  (UNTT=1) 

*  PRINT  *,DEGN 

*  DO  I=0,DEGN 

*  PRINT  *,NUM(I) 
ENDDO 

*  PRINT  *,NO_ FACTS 

*  DO  1=1  ,NO_  FACTS 

*  PRINT  *,DEGD(I) 

*  DO  J=0,DEGD(I) 

*  PRINT  *,DEN(J,  I) 

*  ENDDO 

*  PRINT  *,MULTS(I) 

*  ENDDO 
RETURN 
END 


* 


***************************************************************** 

*  Department  of  Electrical  and  Computer  Engineering      * 

*  Kansas  State  University  * 

*  * 

*  VAX  FORTRAN  source  filename:         SPEC. WRITE. FOR  * 
***************************************************************** 

* 

SUBROUTINE 

SPEC. WRITE  (NUM,  N_DEG,  NO. FACTS,  FACTOR, 

D_DEG,MULTS) 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


ROUTINE: 


DESCRIPTION: 


DOCUMENTATION 
FILES: 


ARGUMENTS: 


NUM 


N  DEG 


NO  FACTS 


FACTOR 


D  DEG 


MULTS 


RETURN: 


This  program  writes  out  the  factored 
form  of  a  rational  function  to  a  file 
specified  by  the  user. 


None. 


The  following  arguments  are  passed  to 
the  subroutine: 

(input)   real 

is  an  array  containing  the  coefficients 

of  the  numerator  polynomial. 

(input)  integer 

is  the  degree  of  NUM 

(input)    integer 

is  the  number  of  factors  in  the 

denominator  polynomial. 

(input)   real 

is  an  array  containing  the 
coefficients  of  each  factor  in 
the  denominator  polynomial. 

(input)    integer 
is  an  array  specifying  the 
degree  of  each  factor  in 
the  denominator  polynomial. 

(input)    integer 
is  an  array  specifying  the 
multiplicity  of  each  factor  in 
the  denominator  polynomial. 

Not  used. 


* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
**************************************************************** 

SUBROUTINE  SPEC. WRITE  ( NUM,  N_ DEG ,  NO.  FACTS,  FACTOR,  D_  DEG ,  MULTS 


ROUTINES 
CALLED: 


AUTHOR: 


DATE  CREATED: 


REVISIONS: 


None. 

James  F.  Stafford 

30Jun88  Version  1.0 

None. 


IMPLICIT 
INTEGER 
REAL*8 
CHARACTER*  15 


NONE 

I,  J,  N_DEGfD_DEG(*)  , NO_ FACTS,  MULTS  (*) 
NUM(0:*)  ,FACTOR(10,0:2) 
FILENAME 
PRINT  *,'  Enter  filename.' 
READ    (*,200)    FILENAME 
200  FORMAT    (A15) 

OFEN    (UNTT=1,   FILE=FILENAME,    STATUS^NEW' ) 
WRITE    (1,*)    N_DEG 
DO  1=0,  N_  DEG 

WRITE    (1,*)    NUM(I) 
ENDDO 

WRITE    (1,*)    NO_ FACTS 
DO  1=1,  NO_  FACTS 

WRITE    (1,*)    D_DEG(I) 
DO  J=0,D_DEG(I) 


WRITE  (1,*)  FACTOR  (I ,  J) 

ENDDO 

WRITE  (1,*)  MULTS(I) 
ENDDO 

CLOSE    (UNTT=lf    STATUS=  *  KEEP' ) 
RETURN 
END 


***************************************************************** 


* 

* 

* 

* 

****** 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


Department  of  Electrical  and  Computer  Engineering      * 
Kansas  State  University  * 

* 

VAX  FORTRAN  source  filename:         TRANSFER. FOR  * 

*********************************************************** 


ROUTINE: 


DESCRIPTION: 


DOCUMENTATION 
FILES: 


ARGUMENTS: 


DEN 


DEGD 


There  are  actually  three  programs 

in  this  file: 

SUBROUTINE 

GETD  ( J,  A,  DEGA,  DEN,  DEGD) 

SUBROUTINE 

GETX  (J,  K,  F,  DEGF,  X,  DEGX) 

SUBROUTINE 

PUTX  (J,  K,  F,  DEGF,  X,  DEGX) 


These  programs    are  a  substitute  for 
a  more  sophisticated  data  structuring 
method.     They  copy  the  coefficients 
for  a  polynomial  embedded  in  a 
higher  dimensional  array  into  a  one- 
dimensional  array,  or  vice  versa. 


None. 


The  following  arguments  are  passed  to 
GETD: 

(input)    integer 

is  a  number  representing  which  factor 

of  the  denominator  is  sought. 

(input)   real 

is  a  two-dimensional  array.     DEN  (I,  J) 
represents  the  coefficient  of  the  Ith 
power  of  x  in  the  Jth  factor  of  the 
denominator  polynomial. 

(input)    integer 

is  an  array.     DEGD(I)   represents  the 
degree  of  the  Ith  factor  in  the 
denominator  polynomial. 

(output)    real 


* 

* 
* 
* 
* 

* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


DEGA 


J,K 


DEGX 


DEGF 


J,K 


is  an  array  to  receive  the  coefficients 
of  the  factor  polynomial. 

(output)    integer 

is  the  degree  of  the  factor  polynomial. 

The  following  arguments  are  passed  to 
GETX: 

(input)  integer 

are  the  coordinates  of  the  term  in  the 

partial  fraction  expansion  that  is 

sought.  See  the  description  of  X. 

(input)   real 

is  a  three-dimensional  array.     X(I, J,  K) 
represents  the  Ith  coefficient  of  the 
numerator  of  the  (J, K)th  term  in  the 
partial  fraction  expansion.     Namely, 
that  term  with  the  Jth  factor  of  DEN 
to  the  Kth  power  as  denominator. 

(input)    integer 

is  an  array.     DEGX(I,J)    represents  the 
degree  of  the  numerator  of  the  (I,  J)th 
term  in  the  partial  fraction  expansion. 
See  the  description  of  X. 

(output)   real 

is  an  array  to  receive  the  coefficients 

of  the  desired  polynomial. 

(output)    integer 

is  the  degree  of  the  polynomial,   F. 

The  following  arguments  are  passed  to 
PUTX: 

(input)    integer 

are  the  coordinates  of  the  term  in  the 
partial  fraction  expansion  that  is  to 
be  updated.   See  the  description  of  X. 

(input)   real 

is  a  three-dimensional  array.     X(I,J,  K) 
represents  the  Kth  coefficient  of  the 
numerator  of  the  (I,  J)th  term  in  the 
partial  fraction  expansion.     Namely, 
that  term  with  the  Ith  factor  of  DEN 
to  the  Jth  power  as  denominator. 


* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 

* 
* 

* 
* 
* 
* 

* 
* 
* 
** 


DEGX 


DEGF 


RETURN: 


ROUTINES 
CALLED: 


AUTHOR: 


DATE  CREATED: 


REVISIONS: 


(input)    integer 

is  an  array.     DEGX(I,J)    represents  the 
degree  of  the  numerator  of  the  (I,J)th 
term  in  the  partial  fraction  expansion. 
See  the  description  of  X. 

(input)   real 

is  an  array  to  containing  the  coefficients 
of  the  polynomial  to  be  embedded  in  the 
partial  fraction  matrix. 

(input)  integer 

is  the  degree  of  the  polynomial,  F. 

Not  used. 


None. 


James  F.  Stafford 


6JunB7     Version  1.0 


None. 


************************************************************** 
SUBROUTINE  GETD(Jf  A,  DEGAf  DEN,  DEGD) 


IMH.ICTT 


NONE 


INTEGER  J,K,DEGA,DEGD(*) 

REAL*8  A(0:*),DEN(0:2,*) 

DEGA=DEGD(J) 
DO  K=0,DEGA 

A(K)=DEN(K,J) 
ENDDO 


RETURN 


END 


SUBROUTINE  GETX(J,K,  F,DEGF,  X,  DEGX) 

IMPLICIT  NCNE 

INTEGER  J,K,L,DEGF,DEGX(10,*) 

REAL*8  F(0:*),X(0:1,10,*) 

DEGF=DEGX(J,K) 

DO  L=0,DEGF 

F(L)=X(LfJ,K) 
ENDDO 
RETURN 
END 

SUBROUTINE 
IMPLICIT 
INTEGER 
REAL*8 

DEGX(J,K)=DEGF 
DO  L=0,DEGF 

X(L,J,K)=F(L) 
ENDDO 
RETURN 
END 

SUBROUTINE 
IMPLICIT 
INTEGER 
REAL*8 


PUTX  ( J,  K,  F,  DEGF,  X,  DEGX) 

NONE 

J,  KfL,  DEGF,  DEGX  ( 10  ,*) 

F(0:*),X(0:1,10,*) 


GET(A,DEGA,B,DEGB) 
NONE 

I,DEGA,DEGB 
A(0:*),B(0:*) 


DEGA=DEGB 
DO   I=OfDEGA 

A(I)=B(i; 
ENDDO 
RETURN 
END 
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Abstract 

The  usual  method  for  solving  linear,  constant-coefficient,  differential  equations 
involves  use  of  the  Laplace  transform.  The  most  difficult  step  in  this  method  of 
solution  is  computing  the  inverse  Laplace  transform  of  a  rational  function.  The 
object  of  this  thesis  is  to  describe  an  algorithm  for  solving  large  systems  of  this 
kind.  The  thesis  demonstrates  that  the  problem  of  solving  such  systems  can  be 
treated  completely  algebraically  once  the  denominator  of  the  rational  function  is 
factored.  It  is  shown  that  the  number  of  operations  required  to  compute  a  suitable 
partial  fraction  expansion  of  a  rational  function  can  be  reduced  by  factoring  the 
denominator  into  irreducible  linear  and  quadratic  factors  in  R[x].  Applications  to 
control  theory  are  discussed.  The  algorithms  are  derived  with  mathematical  rigor. 
Working  FORTRAN  programs  implementing  the  derived  algorithms  are  given  in 
an  appendix.  Electrical  engineers  solving  practical  problems  in  circuit  analysis  and 
control  theory  might  find  these  algorithms  useful. 


